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COMBINED  REFLECTION  AND  DIFFRACTION  BY  A  VERTICAL  WEDGE 


PART  I :  INTRODUCTION 

1.  The  boundary  value  problem  of  linear  wave  reflection  and  diffraction 
by  a  vertical  wedge  of  arbitrary  wedge  angle  has  been  well  formulated  and 
presented  by  Stoker  (1957)  among  many  other  investigators.  The  technique  to 
obtain  an  analytical  solution  for  the  problem  is  also  depicted  in  the  cited 
book.  However,  analytical  solutions  are  not  available  for  the  problem,  except 
for  the  special  case  of  wave  diffraction  by  a  thin  semi-infinite  breakwater, 
that  is,  a  wedge  with  wedge  angle  equal  to  zero. 

2.  The  solution  of  the  thin  semi-infinite  breakwater  was  presented  in 
the  dimensionless  diffraction  diagrams  by  Wiegel  (1962).  The  diagrams  have 
been  especially  useful  in  preliminary  engineering  design  and  have  been  in¬ 
cluded  in  the  Shore  Protection  Manual  (SPM)  (1984).  Although  equally  useful, 
the  combined  reflection  and  diffraction  diagrams  are  not  available,  perhaps 
because  of  the  complexity  of  the  diagrams  which  makes  them  difficult  to  create 
without  using  modern  high-speed  computers  for  computation  and  graphing. 

3.  The  objectives  of  the  present  study  are  (a)  to  obtain  an  analytical 
solution  for  the  combined  wave  reflection  and  diffraction  by  a  vertical  wedge 
of  arbitrary  wedge  angle  subject  to  excitation  of  a  plane  simple  harmonic  wave 
train  coming  from  infinity  and  (b)  to  provide  the  combined  reflection  and  dif¬ 
fraction  diagrams.  The  diagrams  included  in  this  report  have  two  cases:  one 
for  a  thin  semi-infinite  breakwater  and  the  other  for  a  90-deg  vertical  wedge. 
Subroutine  WEDGE  for  computing  the  combined  reflection  and  diffraction  by  a 
vertical  wedge  of  arbitrary  wedge  angle  is  also  documented  in  the  report 
(Appendix  A) . 


PART  II:  BOUNDARY  VALUE  PROBLEM 
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I  Mathematical  Formulation 


4.  In  this  study  our  primary  interest  is  the  wave  reflection  and 
diffraction  by  a  vertical  wedge  of  arbitrary  wedge  angle  in  a  constant  water 
depth  h  *  subject  to  the  excitation  of  monochromatic  incident  waves  of  infi¬ 
nitesimal  amplitude  coming  from  infinity.  Let  (r,9,z)  be  cylindrical  coordi¬ 
nates,  with  z  =  0  representing  the  undisturbed  water  free  surface  and  upward 
direction  representing  the  positive  z-axis.  The  tip  of  the  wedge  is  chosen  to 
be  the  origin  of  the  coordinates  and  two  rigid  walls  of  the  wedge  to  coincide 
with  0=0  and  0=0^,  respectively,  as  illustrated  in  Figure  1.  Carte¬ 
sian  coordinates  (x,y,z),  corresponding  to  the  cylindrical  coordinates,  are 
also  occasionally  used  and  shown  in  the  same  figure.  Therefore,  the  wedge 


angle  Is  -  9^  ,  and  the  water  region  Is  defined  by  0^  ^  6  £  0  and 
0  >  z  >  -h  . 

5.  The  velocity  field  for  the  wave  reflection  and  diffraction  in  an 
ideal  fluid  can  be  represented  by  the  velocity  potential  function  4>(r,0,z,t) 
which  must  satisfy  the  Laplace  equation,  where  t  is  the  temporal  coordinate. 
We  assume  that  the  waves  are  sinusoidal  in  time  with  radian  frequency  to  . 
Water  depth  is  constant,  and  the  bottom  is  rigid  and  impermeable.  Therefore, 
the  vertical  and  temporal  components  of  the  velocity  potential  function,  which 
follow  from  separation  of  variables,  can  be  factored  out  and  the  velocity  po¬ 
tential  written  as 


A  ,  cosh  k(z  +  h)  , „  iwt 

»(r,e,2,t)  -  Ao  - -  ♦<r,9). 


where 

Ao  =  “i8ao/u 
i  « 

g  =  gravitational  acceleration 
aQ  =  incident  wave  amplitude 
k  =  wave  number 

<t>  =  horizontal  component  of  the  velocity  potential  function 

6.  Substituting  Equation  1  into  the  Laplace  equation  and  using  both 
the  kinematic  and  dynamic  boundary  conditions  at  the  free  surface,  the  Laplace 
equation  is  then  reduced  to  the  Helmholtz  equation  which  is  written  in  polar 
coordinates  as  follows: 


2  d2<p  dip  d2<p  ,.22.  n 
r  ^  +  r^  +  ^  +  kr*-° 


where  k  must  be  a  real  number  and  satisfy  the  dispersion  relationship 


u  =  gk  tanh  kh 


7.  The  free  surface  displacement  h  from  the  mean  water  level  z  =  0 
can  be  obtained  from  linear  wave  theory  and  is  represented  as 


8.  Thus  only  the  horizontal  part  of  the  velocity  potential  function  4> 
is  needed  to  be  determined  as  a  solution  of  Equation  2  in  the  water  region 
0^  >  0  >  0  ,  with  the  following  boundary  conditions  at  the  rigid  and  imperme¬ 
able  walls  of  the  wedge: 


34> 

96 


at  6=0  and  6 


(5) 


9.  A  condition  at  infinity  is  also  required  to  ensure  a  unique  solu¬ 
tion.  The  classic  approach  is  to  use  the  Sommerfeld  radiation  condition  at 
infinity  which  states  that  the  scattered  wave  <j>g  must  behave  like  a  cylin¬ 
drical  outgoing  progressing  wave  at  infinity  such  that  ‘ 


,  9r 


+  iktji 


(6) 


The  total  wave  represented  by  <p  is  the  linear  superposition  of  an  incident 


£  ,  a  reflected  wave  from  the  the  9=0  wall  of  the  wedge 


wave 

the  scattered  wave 


and 


4>  from  the  tip  of  the  wedge. 


♦i  +  ^r  +  *s 


(7) 


Equation  6  can  be  satisfied  if 


-ikr 


/kr 


at 


r-*» 


(8) 


10.  The  incident  wave  coming  from  a  large  distance  from  the  tip  of  the 
wedge  is  assumed  to  be  a  plane  progressive  wave  of  amplitude  aQ  and  incident 
angle  a  to  the  x-axis  as  given  by 


ikr  cos  (9-a) 

=  e 


(9) 


Consequently,  the  perfectly  reflected  wave  from  the  y  =  0  wall  of  the  wedge 


Thus  ;  boundary  value  problem  (in  which  the  governing  equation  is  Equa¬ 
tion  2,  the  boundary  condition  is  Eauation  5,  and  the  radiation  condition  is 
Equation  6)  is  completely  formulated. 


Analytical  Solution 

11.  Analytical  solution  to  the  problem  formulated  in  the  preceding 
section  is  obtained  by  following  the  solution  technique  by  Stoker  (1957).  To 
obtain  the  solution,  the  water  region  is  divided  into  three  subregions — I,  II, 
and  III — by  the  incident  wave  ray  passing  through  the  tip  of  the  wedge  and  the 
reflected  wave  ray  reflected  away  from  the  tip  of  the  wedge,  as  shown  in  Fig¬ 
ure  2.  Obviously,  the  total  wave  in  subregion  I  is  the  sum  of  the  incident, 
reflected,  and  scattered  waves;  the  total  wave  in  subregion  II,  where  the 
reflected  wave  does  not  exist,  is  the  sum  of  the  incident  and  scattered  waves; 
and  the  total  wave  in  subregion  III,  where  the  incident  and  reflected  waves 
have  been  shaded  out,  is  only  the  scattered  wave.  For  certain  combinations  of 


/ 


Figure  2.  Three  subregions  and  the  wedge 


the  wedge  angle  and  incident  wave  angle,  subregions  II  and  III  may  not  exist 
at  all.  In  general,  the  solution  function  can  be  written  as 


<t>  =  4>o(r,0)  +  4>s(r ,  0) 


(ID 


where 


(r,0)  =  < 


.  +  4> 

1  r 


7T  -  a  >  6  >  0 

iT  +  a>0>Tr-oi 


(12) 


0  >  0  >  v,  +  a 

o 


The  equation  reveals  that 


is  the  sum  of  the  incident  and  reflected  waves 


p,  and  <fi  and  is  a  known  function.  The  scattered  wave 
i  r 


is  the  only  un¬ 


known  function  to  be  determined  in  the  problem.  Nevertheless,  the  total  wave 


instead  of  the  scattered  wave  A  is  the  desired  solution  to  be  obtained 

s 


in  this  study. 

12.  The  solution  for  the  total  wave  <j>  is  pursued.  The  finite  cosine 
transform  of  <f>  ,  denoted  by  <f>  ,  is  introduced  by  the  formula 


V1T 


<f>(kr  ,n)  = 


/ 


n0 


(kr ,  0)  cos  —  d0 


(13) 


where  n  =  0,1,2,...  are  integers,  and  v  is  related  to  the  wedge  angle  as 
defined  by 


0  —  V7T 

o 


(14) 


Applying  the  finite  cosine  transform  and  using  the  boundary  condition  in  Equa¬ 
tion  5,  Equation  2  becomes 


2  d2<t>  ,  d<p  . 

r  a  2  9r 
9r 


(kr)  ‘ 


-  (av)‘ 


♦  =  0 


(15) 


Equation  15  is  a  form  of  the  Bessel  equation  for  which  general  solutions  are 


the  Bessel  functions  of  the  first  and  second  kinds,  J  ,  (kr)  and  Y  ,  (kr)  , 

n/ v  n/ v 


respectively.  Since  Y^^(kr)  are  singular  at  the  origin,  the  solution  is 


chosen  to  be 


where  a  are  constants  to  be  determined, 
n 

13.  Taking  the  finite  cosine  transform  of  Equation  11  and  using  Equa¬ 
tion  16,  we  have 


d>  cos  —  d0  =  a  J  ,  (kr) 
s  v  n  n/v 


V1T 


n0 

4  cos  —  d0 
Yo  v 


(17 


or 


s 


a 

n. 


n/v 


(kr)  -  4o 


(18 


Then  applying  the  operation  lim^^  /r(9/9r  +  ik)  to  both  sides  of  Equa¬ 
tion  18,  and  using  the  Sommerfeld  radiation  condition  (Equation  6)  we  have 


lim 

r-K» 


/^l 


(fe  *  “) 


a„  J„/v<kr) 


n0 

d>  cos  —  d0 

T0  V 


0 


(19 


14.  Equation  19  can  be  asymptotically  evaluated  to  determine  a^  . 

Firstly,  the  first  term  involving  the  Bessel  function  is  evaluated.  The  func 

tion  J  .  (kr)  at  r-*»  behaves  asymptotically  (Abramowitz  and  Stegun  1964) 
n  /  v 

as  follows: 


Jn/v(kr) 


(20 


Hence,  we  have 


lim  /?/!-  +  ik\ J  ,  (kr) 


2k  i(kr-nTr/2v-Hr/4) 
\  —  e 


(21 


I 

I 

■ 

Secondly,  the  second  term  involving  the  integral  of  is  evaluated.  The 

asymptotic  behavior  of  the  integral  over  9  =  (O.vtt)  and  at  large  distance 
I  r can  be  found  by  the  method  of  stationary  phase.  The  integral,  after  sub¬ 
stituting  from  Equations  9,  10,  and  12,  can  be  written  as 


V7T  TT-Ct 

J  4>o  cos  ^  d6  =y  [ 


ikr  cos(9-a)  ikr  cos(0+a)l  n9 
e  +  e  J  cos  — 


d6 


n+a 


I 

Tr-ot 


.  ikr  cos(9-a)  n9  R 
+  /  e  cos  —  d9  (22) 


In  the  integrals,  there  are  three  points  of  stationary  phase  at  0  =  a  and 
0  =  it  ±  a  .  If  the  same  argument  as  that  of  Stoker  (1957)  is  followed,  of  the 
three  contributions  only  the  first  one  0  =  a  furnishes  a  nonvanishing  con¬ 
tribution  for  r-*°°  when  the  operator  >^r(9/3r  +  ik)  is  applied  to  it.  The 
physical  significance  of  this  statement  is  that  only  the  incident  wave  is 
effective  in  determining  the  cosine  coefficients  of  the  solution.  Therefore, 


(>  cos  ~  d0  ~  2/2iTk  cos  —  e 
o  v  v 


na  i(kr+ir/4) 


(23) 


Substituting  Equations  21  and  23  into  Equation  19,  we  obtain  the  unknown  co¬ 
efficients  a  : 

n 


a 

n 


na  inir/2v 
2tt  cos  —  e 


(24) 


15.  Since  the  solution  <t>  in  the  cosine  series  expression  is 


the  solution  is  obtained  by  substituting  Equations  16  and  24  into  Equation  25 
as  follows: 


<t>(r ,  0)  =  ^ 


JQ(kr) 


+  2 


I 


inn/2vr  ,  na  n0 

e  J  ,  (kr)  cos  —  cos  — 

n/v  v  v 


n=l 


(26) 


Equation  26  is  the  solution  for  the  combined  wave  reflection  and  diffraction 
by  a  vertical  wedge  of  arbitrary  wedge  angle  and  is  considered  to  be  extended 
from  the  solution  by  Stoker  (1957)  who  only  solved  the  problem  of  a  thin  semi¬ 
infinite  breakwater.  The  solution  in  Equation  26  and  the  one  by  Stoker  are 
not  only  in  nonclosed  form  but  also  in  terms  of  Bessel  functions.  It  seems 
that  the  calculations  of  the  solutions  are  very  difficult  without  using  a  mod¬ 
ern  high-speed  computer.  This  is  probably  the  reason  why  Stoker  arrived  at 
his  solution  expressed  in  the  same  cosine  series  but  did  not  use  it  to  calcu¬ 
late  the  result.  Instead,  he  further  transformed  the  expression  into  a  very 
complex  integral  form  for  further  approximation  in  calculating  the  result. 

16.  Notably,  the  solution  at  the  origin  point  is  obtained  by  simply 
substituting  r  *  0  into  Equation  26  to  arrive  at 

♦(0,e>-|  (27) 

Therefore,  wave  response  at  the  origin  point  depends  only  on  the  wedge  angle 
and  does  not  depend  on  the  incident  wave  angle. 

Two  Special  Cases 

17.  The  solutions  for  two  special  cases  are  used  to  verify  Equation  26: 
one  for  the  case  of  a  thin  semi-infinite  breakwater  and  the  other  for  the  case 
of  an  infinite  wall  extending  from  x  =  -®  to  00  . 

18.  The  vertical  wedge  should  reduce  to  a  thin  semi-infinite  breakwater 

as  the  wedge  angle  reduces  to  0  deg.  Therefore,  solution  of  the  combined  wave 
reflection  and  diffraction  by  a  thin  semi-infinite  breakwater  is  obtained  by 
substituting  v  =  2  (that  is,  0  ■  2tt)  into  Equation  26  which  then  becomes 


<Kr,6)  =  J0(kr)  +  2  £ 

n=l 


J  /0(kr)  cos  ■=-  cos  zr- 
n  1 1  i  i 


Equation  28  is  precisely  the  same  one  obtained  by  Stoker  (1957). 

19.  The  vertical  wedge  should  also  become  an  infinite  wall  extending 
from  x  =  -00  to  00  with  the  water  occupying  only  the  half  plane  of  y  2  0  as 
the  wedge  angle  increases  to  180  deg.  In  this  situation  the  scattered  wave  is 
absent  from  the  solution,  and  the  total  wave  is  only  the  sum  of  the  incident 
and  reflected  waves  as  follows: 

<j>(r,6)  =  eikr  co^e~a>  +  eikr  c°s(0-Kx)  (29) 

After  expansion  of  the  exponential  functions  in  terms  of  Bessel  functions 
(Abramowitz  and  Stegun  1964) ,  Equation  29  becomes 


,0)  =  2  Jg(kr)  +  2  ^  inJn(kr)  cos  not  cos 


Equation  30  is  the  same  equation  reduced  from  Equation  26  by  substituting 
v  =  1  into  it. 


PART  III:  CALCULATION  AND  RESULTS 


20.  Results  of  the  combined  reflection  and  diffraction  by  a  wedge  of 
arbitrary  wedge  angle  can  be  calculated  from  Equation  26.  Since  the  solution 
is  not  only  in  terms  of  Bessel  functions  but  also  in  a  nonclosed  form,  the 
computer  program  WEDGE  is  therefore  written  to  calculate  the  solution. 

21.  In  the  program  the  subroutine  BESJ  for  calculating  Bessel  function 
of  fractional  or  integer  order  was  used.  The  subroutine  was  originally  writ¬ 
ten  by  Amos,  Daniel,  and  Weston  in  1975  (Morris  1984)  and  is  collected  in  the 
Naval  Surface  Weapons  Center  Library  of  Mathematics  Subroutines  (Morris  1984). 

22.  In  the  calculation  the  summation  of  the  infinite  terms  in  Equa¬ 
tion  26  was  carried  out  to  the  term  which  is  preceded  by  eight  successive 

terms  of  the  absolute  value  of  the  Bessel  function,  all  equal  to  or  less  than 

—8  —8 
10  .  The  solution  has  a  truncation  error  less  than  10  ,  and  it  is  of  the 

order  of  one. 

23.  In  this  study,  results  of  the  combined  wave  reflection  and  diffrac¬ 
tion  for  the  wedge  are  calculated  for  two  cases:  one  for  a  vertical  wedge  of 
0-deg  wedge  angle  and  the  other  for  a  vertical  wedge  of  90-deg  wedge  angle. 

Vertical  Wedge  of  0-Deg  Wedge  Angle 

24.  When  the  wedge  angle  is  equal  to  zero,  the  wedge  is  actually  a  thin 
semi-infinite  breakwater  extending  from  x  =  0  to  00  .  Figure  3  shows  the 
thin  semi-infinite  breakwater  along  with  the  polar  coordinates.  In  this  case 
the  diffraction  results  for  various  incident  wave  angles  in  the  water  region 
from  0  =  it  to  2t r  and  r/A  <  10  ,  where  A  is  the  incident  wave  length, 
have  already  been  presented  by  Wiegel  (1962)  and  are  shown  in  the  SPM,  Vol¬ 
ume  I  (1984).  The  present  results  combine  reflection  and  diffraction  effects 
and  cover  the  water  region  from  9  =  0  to  2tt  and  r/A  <  10  .  Therefore, 
the  present  results  for  this  particular  case  can  be  considered  to  be  a  comple¬ 
mentary  and  extended  version  to  the  ones  in  the  SPM. 

25.  In  this  study  wave  response  was  calculated  at  1,460  grid  points 

intersected  by  r/A  =  0.5, (0.5), 10  ,  which  means  that  the  values  of  r/A  are 
from  0.5  to  10.0  with  each  value  increment  being  0.5.  Hereafter,  all  similar 
expressions  are  to  be  interpreted  in  the  same  way  (e.g.,  0  =  0,  (r/36)  ,2tt  for 


Figure  3.  Thin  semi-infinite  breakwater  and  polar  coordinates 

the  incident  wave  angle  a  -  0,(tr/12)  ,ir)  .  The  wave  response  at  the  origin 
point  is  obtained  substituting  v  =  2  into  Equation  27,  as  follows: 

♦(0,0)  =  1  (31) 

Those  calculated  values  were  used  to  interpret  the  value  for  each  non¬ 
overlapping  pixel  of  size  O.lr/A  by  O.lr/A  in  the  area  within  the  lOr/A 
radius  from  the  origin.  A  diagram  was  then  constructed  by  patching  those 
pixels  over  the  entire  area.  The  wave  response  diagrams  for  each  incident 
wave  angle  are  shown  in  Figures  4  through  15.  Notably,  the  values  in  the  dia¬ 
grams  constitute  the  amplification  factor  which  is  defined  as  the  ratio  of  the 
total  wave  height  to  the  incident  wave  height.  Therefore,  in  subregions  II 


Figure  4.  Amplification  factor  diagram  for  the  thin  semi-infinite 
breakwater  for  incident  wave  angle  **  0  deg 


and  III  (as  defined  in  Figure  2)  where  the  reflected  wave  is  absent,  the  am¬ 
plification  factor  is  essentially  the  diffraction  coefficient  as  defined  in 


the  SPM. 


26.  Figures  4  through  15  reveal  that  the  amplification  factors  in  sub- 
region  I  change  very  rapidly  between  0  and  2.35  over  the  subregion,  and  the 
diagram  patterns  become  very  complex  because  of  the  interesting  superposition 
of  the  incident,  reflected,  and  scattered  waves.  (In  the  legend  of  Figures  4 
through  15,  the  width  of  the  pixel  is  one  incident  wave  length,  and  the  values 
are  amplification  factors.)  Such  patterns  would  be  very  difficult  to  con¬ 
struct  without  using  a  high-speed  computer  and  computer  graphics.  In  sub- 
regions  II  and  III,  the  amplification  factors  change  smoothly  from  1.15 
roughly  along  the  reflected  wave  ray  reflected  from  the  origin  point  to  nearly 
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Figure  13.  Amplification  factor  diagram  for  the  thin  semi-infinite 
breakwater  for  incident  wave  angle  =  135  deg 
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0.00  at  the  back  wall  of  the  wedge  in  the  shadow  zone.  The  diagram  patterns 
are  relatively  smooth  and  simple. 

27.  Notably,  the  diagrams  do  not  include  phase  information  of  the  wave 
response  which  is  usually  unimportant  in  most  engineering  practice.  Should 
the  phase  of  the  wave  response  need  to  be  known,  one  can  always  use  the  com¬ 
puter  program  WEDGE  to  calculate  it. 

28.  The  contour  diagrams  of  the  amplification  factor,  similar  to  the 
ones  presented  by  Wiegel  (1962)  and  shown  in  the  SPM  (1984),  were  also  plotted 
as  typically  shown  in  Figure  16.  Examination  of  those  contour  diagrams  indi¬ 
cates  that,  in  subregion  III,  the  results  are  identical  to  Wiegel' s  results. 
But  in  subregion  II  the  contour  patterns  for  the  amplification  factor  (or  the 
diffraction  coefficient  K'  used  in  the  SPM  (1984))  equal  to  1.0;  thus  the 
present  results  are  far  more  complicated  than  Wiegel' s.  The  author  believes 
that  Wiegel 's  results  may  lose  accuracy  because  of  insufficient  resolution  of 
the  computational  tools  during  the  late  fifties  and  early  sixties.  Neverthe¬ 
less,  such  inaccuracies  are  usually  either  tolerable  or  immaterial  in  most 
engineering  practice. 

29.  The  contour  patterns  in  subregion  I  are  very  complex,  and  it  is 
difficult  to  track  specific  contours.  Therefore,  for  clarity  only  the  patched 
diagrams  are  presented,  and  the  contour  diagrams  are  omitted  in  this  report. 

Vertical  Wedge  of  90-Deg  Wedge  Angle 


30.  When  the  wedge  angle  is  equal  to  tt/2 ( Qq  =  3n/2)  ,  the  vertical 
wedge  occupies  the  entire  fourth  quadrant  of  the  space  as  shown  in  Figure  17. 
Wave  response  was  calculated  at  1,100  grid  points  intersected  at  r f\ 
m  0.5,  (0.5)  ,  10.0  and  0  ■  0,  (ir/36)  ,3tt/2  for  the  incident  wave  angle  a 
«  0,(w/12),tt  .  The  wave  response  at  the  origin  is  obtained  by  substituting 
v  *  1.5  into  Equation  27  as  follows: 

♦(0,9)  -  \  (32) 

Those  calculated  results  were  used  to  construct  the  amplification  factor  dia¬ 
grams  by  following  the  same  procedures  described  for  the  case  of  the  thin 
semi-infinite  breakwater.  The  diagrams  are  shown  in  Figures  18  through  27. 
Because  of  symmetry  of  the  results,  the  diagrams  for  the  incident  wave  angle 
a  >  3it/4  can  be  obtained  from  those  for  an  incident  wave  angle  of  tt  -  a  . 


Figure  16.  Amplification  factor  contour  diagram  for  the  thin 
semi-infinite  breakwater  for  incident  wave  angle  =  90  deg 


Therefore,  the  diagrams  for  the  incident  wave  angles  a  *  tt/6 ,  (it/12)  , it  are 
omitted  in  this  report. 

31.  The  diagrams  indicate  that,  for  each  corresponding  incident  wave 
angle,  the  results  in  subregion  I  are  very  similar  to  those  obtained  for  the 
vertical  wedge  of  0-deg  wedge  angle.  However,  the  results  in  subregions  II 
and  III  from  both  wedges  are  discernibly  different. 


O-  tot.' 


24.  Amplification  factor  diagram  for  the 
g  wedge  for  incident  wave  angle  =  90  deg 
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PART  IV:  CONCLUSION 


32.  An  analytical  solution  for  the  combined  wave  reflection  and  dif¬ 
fraction  by  a  vertical  wedge  of  arbitrary  wedge  angle  is  obtained  and  ex¬ 
pressed  in  Equation  26.  The  analytical  solution  is  in  terms  of  Bessel 
functions  and  in  nonclosed  form.  The  computer  subroutine  WEDGE,  written  for 
calculating  the  solution,  is  documented  in  Appendix  A. 

33.  The  amplification  factor  diagrams  for  a  vertical  wedge  of  0-deg 
wedge  angle  and  a  vertical  wedge  of  90-deg  wedge  angle  are  calculated  and  pre¬ 
sented.  The  calculated  results  indicate  that  the  wave  response  in  subre¬ 
gion  I,  where  the  incident,  reflected,  and  scattered  waves  all  exist,  is  in  a 
very  complex  pattern  with  the  amplification  factor  varying  from  2.35  to  0.0 
over  the  subregion.  The  wave  response  in  subregions  II  and  III  is  a  rela¬ 
tively  simple  pattern  with  the  amplification  factor  decreasing  from  1.15 
roughly  along  the  reflected  wave  ray  reflected  from  the  origin  point  to  nearly 
0.00  at  the  back  wall  of  the  wedge  in  the  shadow  zone. 

34.  Diagrams  of  the  special  case  of  a  vertical  wedge  of  0-deg  wedge 
angle  can  be  considered  complementary  and  extended  versions  to  the  ones  pre¬ 
sented  in  the  SPM  (1984). 
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APPENDIX  A:  SUBROUTINE  WEDGE 


1.  Subroutine  WEDGE  is  used  to  calculate  the  value  of  <t>  in  Equa¬ 
tion  26  which  is  generally  a  complex  number.  Its  absolute  value  is  the  ampli¬ 
fication  factor,  and  its  phase  is  the  phase  indicator  from  the  phase  of  the 
incident  wave.  As  mentioned  in  the  main  text,  <f>  is  a  function  of  the  Bessel 
function  of  either  fractional  or  integer  order,  depending  on  the  wedge  angle, 
and  is  the  summation  of  a  series  of  infinite  terms.  The  subroutine  BESJ,  doc¬ 
umented  in  the  Naval  Surface  Weapons  Center  (NSWC)  Library  of  Mathematics  Sub¬ 
routines  (Morris  1984)*  is  used  in  the  WEDGE  subroutine.  The  programming  of 
the  WEDGE  subroutine  is  very  straightforward  if  a  truncation  term  in  the  se¬ 
ries  in  Equation  26  is  determined.  The  program  is  written  in  FORTRAN  language 
and  is  listed  in  this  appendix. 


Description 


2.  The  following  subroutine  is  available  for  computing  <f>  in 
Equation  26: 


CALL  WEDGE (P , FABS , FPHA , XRL , XTH , WEDGEA , WAVEA , IDX) 

where  the  arguments  are  all  real  values  except  F  which  is  a  complex  value. 

Input  arguments  are  as  follows: 

a.  (XRL,XTH)*(r/A,0)  where  (r,8)  are  polar  coordinates  of  the  lo¬ 
cation  where  $  is  to  be  computed,  and  A  is  the  incident  wave 
length.  Therefore,  XRL  is  the  radius  vector  or  radius  distance 
normalized  by  the  incident  wave  length.  XTH  is  the  vectorial 
angle  in  degree. 

b.  WEDGEA  =  wedge  angle  in  degree. 

c.  WAVEA  =  incident  wave  angle  in  degree. 

d.  IDX  «  an  index  (set  to  0  in  this  subroutine). 

Output  arguments  are  as  follows: 


*  References  cited  in  the  Appendix  can  be  found  in  the  References  at  the  end 
of  the  main  text. 


b. 

c . 


F  =  <p  in  Equation  26  (wave  response  normalized  by  incident 
wave  amplitude) . 

FABS  =  amplification  factor,  the  absolute  value  of  <p  . 

FPHA  =  phase  difference,  the  phase  of  tj>  . 


Example  and  Test  Run 

3.  To  serve  as  an  example  as  well  as  to  ensure  the  subroutine  is  the 
correct  one,  the  user  should  run  the  test  program  listed  in  Figure  A1  and  make 
sure  the  output  is  the  same  as  that  listed  in  Table  A1 . 


PROGRAM  WEDGE  1 (INPUT, OUTPUT, TAPE5=INF'UT,TAPE6  =  0UT  PUT) 
COMPLEX  F 
PI  =3.  141592654 
WR I  TE  (  6 , 4 ) 


TEST  1 
TEST2 
TEST3 
TEST4 


l  FORMAT (//3X ,  WEDG-ANG  WAV-ANG  LOCATION 

1  ABS-VAL  PHASE ' / 3  X , '  (DEG)  (DEG)  XRL 

AMP-FAC  (RAD)  '  / 3X , 35 ( '  -  '  ) / ) 


2ALIZED) 

WEDGEA=90. 

W  A  V  E  A  =  1 3  5 . 

X  R  L  =  2 . 0 
X  T  H = 3  0 . 

I D  X.  =  0 

CALL  WE DGE(F, FABS, FPHA, XRL, XTH, WEDGE A , WAVEA , IDX ) 
WRITE (6,40)  WEDGEA, WAVEA, XRL, XTH,F, FABS, FPHA 
40  FORMAT (IX, 8F9 . 2) 

STOP 

END 

Figure  A1 .  Computer  program  list  1 


WAVE  RESPONSE  TESTS 
XTH ( DEG )  ( N0RMTEST6 

TEST7 


TEST8 
TEST? 
TEST  1 0 
TEST  1 1 
TEST  1 2 
TEST  1 3 
TEST  1 4 
TESTIS 
TEST  1 6 
TEST  1 7 


Table  A1 

Sample  Output  of  the  Test  Program 


WEDG-ANG  WAV-ANG  LOCATION  WAVE  RESPONSE  ABS-VAL  PHASE 

(DEG)  (DEG)  XRL  XTH (DEG)  (NORMALIZED)  AMP-FAC  (RAD)  ! 

| 


The  input  arguments  in  the  test  programs  are  as  follows: 

a.  WEDGEA  =  90;  the  wedge  angle  is  90  deg. 

b.  WAVEA  =  135;  the  incident  wave  angle  is  135  deg. 


XRL  =  2.0  and  XTH  =  30;  the  location  of  the  wave  response  <f>  to 
be  computed  is  at  radial  vector  of  two  incident  wave  length  dis¬ 
tances  and  vectorial  an6le  of  30  ueg  in  polai.  coordinates. 


d.  IDX  =  0;  the  index  IDX  is  set  to  0. 


This  case  is  shown  in  Figure  A2.  The  outputs  are  given  in  Table  A1  which  is 
self-explanatory . 

5.  If  the  location  of  the  wave  response  <f>  to  be  computed  is  a  very 
large  distance,  for  example  XRL  greater  than  18  or  so,  the  output  might  print 
the  message  that  the  number  of  terms  is  insufficient  for  summation  in  comput¬ 
ing  <f>  .  In  this  situation,  the  user  must  replace  the  integer  200  in 

PARAMETER (NN= 200) ,  listed  in  Card  WEDGE15  in  the  SUBROUTINE  WEDGE  list,  by  a 
larger  integer  to  ensure  the  accuracy. 
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Figure  A2.  Example  problem 
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6.  Subroutine  WEDGE,  which  is  listed  in  this  section,  calls  subroutine 


BESJ  which,  in  turn,  calls  subroutine  JAIRY  and  function  GAMLN.  Subroutines 


PFSJ  and  JATRY  and  function  GAMLN  are  borrowed  from  the  NSWC  Library  of  Mathe¬ 
matics  Subroutines  (Morris  1984).  Including  these  borrowed  subroutines  here 
in  the  list  is  only  for  the  purpose  of  allowing  subroutine  WEDGE  to  be  self- 
contained  and  complete.  The  computer  program  of  subroutine  WEDGE  is  listed  in 
Figure  A3. 
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SUBROUTINE  WEDGE i F , F ABS , FPHA , XRL , X TH , WEDGE A , WAVE A , I DX )  WEGDE1 


c 

* 

*************** 

******** 

* 

*  *  * 

***** 

*  *NOT  I CE 1 

c 

* 

THIS  COMPUTER  PROGRAM  WAS 

WRITTEN  BY  H.S. 

CHEN  OF 

CERC  IN 

♦N0TICE2 

r 

* 

1 936  UNDER  CORF'  CIVIL  WORK  ELD  PROGRhH. 

NEITHER 

ANY  OF 

*NOT 1CE3 

c 

* 

AGENCIES  NOR  ANY  INDIVIDUAL  ASSUMES  ANY  LEGAL  LIABILITY  OR 

♦N0TICE4 

c 

♦ 

RESPONSIBILITY  FOR  THE  ACCURACY  OF  THE  PROGRAM. 

*NOT I CE5 

c 

* 

*************** 

******** 

* 

*  *  * 

***** 

*  *NQT ICE6 

c 

- 

- 

-  -  - 

-  -WEDGE2 

c 

WAVE  REFLECTION  AND  DIFFRACTION  BY  A  VERTICAL 

WEDGE. 

WEDGE3 

c 

WEDGE4 

c 

INPUT  :  XRL  =  R/L,  (RADIUS  VECTOR) / (WAVE 

LENGTH) 

WEDGES 

c 

X  TH  =  VECTORIAL 

ANGLE  IN  DEGREE 

WEDGE6 

r 

u 

WEDGEA  =  WEDGE  ANGLE  IN  DEGREE 

WEDGE7 

C 

WAVEA  =  INCIDENT 

WAVE  ANGLE  IN  DEGREE 

WEDGES 

c 

OUTPUT  :  F  =  PHI,  WAVE 

RESPONSE  NORMALIZED  BY 

INCIDENT 

WEDGE9 

c 

WAVE  AMPLITUDE 

WEDGE10 

c 

FABS  =  ABSOLUTE 

VALUE  OF  PHI,  THE 

AMPLIFICATION 

WEDGE11 

c 

FACTOR 

WEDGE12 

c 

FPHA  =  PHASE  OF 

PHI  IN  RADIAN. 

WEDGE13 

c 

- 

- 

-  -  - 

-  -WEDGE14 

PARAMETER (NN=200) 

WEDGE15 

DIMENSION  BJ (NN) ,W (NN) , XM (NN> 

WEDGE16 

COMPLEX  F , TM 

WEDGE17 

DATA  TOLR/i.E-8/, ITER/8/ 

WEDGE18 

P I =  3 . 141592654 

WEDGE  1 9 

CPI=P I / 1 80 . 

WEDGE20 

XKR=XRL*2 . 0*P  I 

WEDGE21 

TH=XTH*CPI 

WEDGE22 

WA=WAVEA*CP I 

WEDGE23 

XNU= (360. -WEDGEA) / 180. 

WEDGE24 

XM ( 1 ) =  1 . 0 

WEDGE25 

IF(IDX.NE.O)  GOTO  14 

WED6E26 

CALL  BESJ ( XKR , 0 . 0 , 1 , W , NZ ) 

WEDGE27 

BJ ( 1 ) =W ( 1 ) 

WEDGE28 

I C0UNT=0 

WEDGE29 

DO  10  N= 1 , NN 

WEDGE30 

N 1 =N+ 1 

WEDGE31 

I  COUNT  = I  COUNT  + 1 

WEDGE32 

IF ( ICOUNT. LE. ITER)  GOTO  B 

WEDGE33 

NNN  =  N 

WEDGE34 

GOTO  14 

WEDGE35 

8  XM(N1)=FL0AT  <  N ) / X  N  U 

WEDGE36 

M= INT ( XM (N 1 ) ) 

WEDGE37 

ALPHA=XM(N1)-M 

WEDGE38 

M1=M+ 1 

WEDGE39 

CALL  BESJ (XKR, ALPHA, Ml ,W,NZ) 

WEDGE40 

IF(Nl.EQ.NN)  WR I TE ( 6 , 9 )  XKR , ALPHA , Ml , W ( M 1 ), NZ 

WEDGE41 

9  FORHAT ( / '  ****  NO.  OF  TERMS  FOR  SUMMATION  IS 

INSUFFICIENT  **** 

'  ,  WEDGE42 

1  i '  XKR, ALPHA, Ml , W ( M 1 ) ,NZ  ='  , 

2F10.4, I5,E15.6, 

I 

5/) 

WEDGE43 

B J  ( N 1 ) =W ( M 1 ) 

WEDGE44 
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«*. 


10 

14 


30 


IF ( ABS ( BJ (Nil ) .6T.T0LR)  IC0UNT=0 
CONTINUE 
CONTINUE 
F  =  BJ  <  1 ) / 2 . 

DO  20  N= 1 , NNN 
N1=N+ 1 
XMN=XM ( N 1 ) 

TM=<0.0, 1.0)**XMN*BJ (N1)*C0S(XMN*WA)*C0S(XMN*TH) 

F=F+TM 

CONTINUE 

F=4. /XNU*F 

FR=REAL (F) 

FI=AIMAG(F) 

FABS=SQRT (FR*FR+FI*FI) 

IFINA.LE. l.E-B)  FABS=FABS/2. 

IFtFABS.LT. TOLR)  GOTO  30 
FPHA=ATAN2 (FI ,FR) 

RETURN 
FPHA=0 . 0 
RETURN 
END 

SUBROUTINE  BESJ ( X , ALPHA , N , Y , NZ > 

WRITTEN  BY  D.E.  AMOS,  S.L.  DANIEL  AND  M.K.  WESTON, 
REFERENCE  SAND-75-0147 
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ABSTRACT 

BESJ  COMPUTES  AN  N  MEMBER  SEQUENCE  OF  J  BESSEL  FUNCTIONS 
J/SUB ( ALPHA+K- 1 ) / ( X ) ,  K=l,...,N  FOR  NON-NEGATIVE  ALPHA  AND  X 
A  COMBINATION  OF  THE  POWER  SERIES,  THE  ASYMPTOTIC  EXPANSION 
FOR  X  TO  INFINITY  AND  THE  UNIFORM  ASYMPTOTIC  EXPANSION  FOR 
NU  TO  INFINITY  ARE  APPLIED  OVER  SUBDIVISIONS  OF  THE  <NH,X) 

PLANE.  FOR  VALUES  OF  <NU,X)  NOT  COVERED  BY  ONE  OF  THESu 
FORMULAE,  THE  ORDER  IS  INCREMENTED  OR  DECREMENTED  BY  INTEGER 
VALUES  INTO  A  REGION  WHERE  ONE  OF  THE  FORMULAE  APPLY.  BACKWARDBESJ 
RECURSION  IS  APPLIED  TO  REDUCE  ORDERS  BY  INTEGER  VALUES  EXCEPTBESJ 
WHERE  THE  ENTIRE  SEQUENCE  LIES  IN  THE  OSCILLATORY  REGION.  IN  BESJ 
THIS  CASE  FORWARD  RECURSION  IS  STABLE  AND  VALUES  FROM  THE  BESJ 
ASYMPTOTIC  EXPANSION  FOR  X  TO  INFINITY  START  THE  RECURSION  BESJ 
WHEN  IT  IS  EFFICIENT  TO  DO  SO.  LEADING  TERMS  OF  THE  SERIES  ANDBESJ 
UNIFORM  EXPANSION  ARE  TESTED  FOR  UNDERFLOW.  IF  A  SEQUENCE  IS  BESJ 
REQUESTED  AND  THE  LAST  MEMBER  WOULD  UNDERFLOW,  THE  RESULT  IS  BESJ 
SET  TO  ZERO  AND  THE  NEXT  LOWER  ORDER  TRIED,  ETC.,  UNTIL  A  BESJ 
MEMBER  COMES  ON  SCALE  OR  ALL  MEMBERS  ARE  SET  TO  ZERO.  OVERFLOWBESJ 
CANNOT  OCCUR.  BESJ1  CALLS  SUBROUTINE  JAIRY  AND  FUNCTION  GAMLN . BESJ 

BESJ 

DESCRIPTION  OF  ARGUMENTS  BESJ 

BESJ 

INPUT  BESJ 

X  -  X.GE.O  BESJ 

ALPHA  -  ORDER  OF  FIRST  MEMBER  OF  THE  SEQUENCE.  ALPHA. GE.O  BESJ 
N  -  NUMBER  OF  MEMBERS  IN  THE  SEQUENCE,  N.GE.l  BESJ 

OUTPUT  BESJ 

Y  -  A  VECTOR  WHOSE  FIRST  N  COMPONENTS  CONTAIN  BESJ 

VALUES  FOR  J /SUB (ALPHA+K- 1 )/( X ) ,  K=l . N  BESJ 
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-  ERROR  INDICATOR 


NZ  =0 
N  Z  =  —  1 
NZ  =  -2 
N  Z  =  —  3 
NZ.GT.O 


NORMAL  RETURN  -  COMPUTATION  COMPLETED 
X  IS  LESS  THAN  0.0 
ALPHA  IS  LESS  THAN  0.0 
N  IS  LESS  THAN  1 

LAST  NZ  COMPONENTS  OF  Y  SET  TO  0.0 
BECAUSE  OF  UNDERFLOW 


ERROR  CONDITIONS 


IMPROPER  INPUT  ARGUMENTS  -  A  FATAL  ERROR 
UNDERFLOW  -  A  NON-FATAL  ERROR  (NZ.GT.O) 


DOUBLE  PRECISION  DX , TRX , DTM , DFN 
DIMENSION  Y (N) 

DIMENSION  C (11,10) .ALFA (26, 4) .BETA (26, 5) 

DIMENSION  Cl (88) , C2 ( 22 ) 

DIMENSION  A l  (52)  ,A2(52)  ,B1 (52) , B2 (52 !  , B3 (26 ) 

DIMENSION  GAMA (26) , TEMP (3) , KMAX (5 ) , AR (8 ) , BR ( 10) , UPOL (10) 
DIMENSION  FNULIM (2) ,PP<4) 

DIMENSION  CR(10)  , DR ( 1 0 ) 


EQUIVALENCE  ( C ( 1 , 1 ) , C 1  ( 1 )  ) 
EQUIVALENCE  (C < 1 ,9) , C2 ( 1 > > 
EQUIVALENCE  (ALFA ( 1 , 1 ) , A1 ( 1 ) ) 
EQUIVALENCE  ( ALFA ( 1 , 3 )  ,  A2  ( 1 ) ) 
EQUIVALENCE  (BETA  ( 1 , 1 ) , B 1  ( 1 ) ) 
EQUIVALENCE  ( BETA ( 1 , 3 )  ,  B2  ( 1 ) ) 
EQUIVALENCE  ( BETA ( 1 , 5 ) , B3  ( 1 ) ) 

DATA  EL  I M 1 , tL I  M2 , TOL  /  667. 


l.E-15 


DATA  PP(l)/8. 72909 15393555E  +  00/ ,  PP  (2) il.  6569393226503E-01 / , 
1  PP (3 ) /  1 . 2457857686559E— 01 / ,  PP<4)/7. 70 13374743 039E- 04/ 


TQLS=LN  < 1 . E — 3 ) 
DATA  TOLS 


/-6.9077552789821E+00/ 


DATA  UPOL(l) ,C0N1 ,C0N2,CQN3,C0N548  /  1 . OOOOOOOOOOOOOE+OO , 

1  6.6666666666667E-01,  3.3333333333333E-01  ,  1 . 4 1 42 1 3562373 1 E+00 , 

2  1.0416666666667E-01/ 

DATA  RTWO,PDF,RTTP,PIDT  /  1 . 3483997249265E  +  00 , 

1  7. 85398 1633974 5E -01,  7 . 9788456080286E-0 1 ,  1 . 5707963267949E+00/ 

DATA  FNULIM ( 1 ) / 100. / ,  FNULIM (2) /60. / 

CE=-ALQG ( TOL )  ,  TCE=-0. 75*AL0G (TOL) 

DATA  CE  ,  TCE  /  3. 453877639491 1E+01 ,  2. 59040822961 83E+01  / 

DATA  INLIM  /  150  / 

DATA  AR (  1 )  /  8.3550347222222E-02/,  AR ( 2 ) / 1 . 2B22657455633E-0 1 / , 
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AR(3i / 2. 9  1349026464  14E-01/.  AR ( 4 ) /8 . 8 1 62726 744376E-0 1 / ,  BESJ189B 

A  R  (  5  >  /  3 . 3  2  1  4  0  8  2  8 1  B  6  2  8  E  +  0  0 / ,  AR  (6> /  1 . 4995762986863E+01 / ,  BE  S J 1 89  C 
AR (7) /?. 892301301  15B6E+01/ ,  AR < 8 ! / 4 . 7445 1 53886826E+02/  BESJ 169D 

BESJ190 

DATA  BRU)  /  -  1 .  4  5  8  3  3  3  3  3  3  3  3  3  3  E  -  0  1  /  ,  BR(2)  /-9. 8741319444444E-02/ ,  BSE  J  1 9 1 A 

BR  (3)  - 1 . 4  3  3 1 2  0  5  3  9 1 5  9  0  E  -  0 1  /  ,  BR  ( 4)  / -3.  172272026784  IE-01  /  ,  BSE  J 1 9 1 B 

BR  (5 )  / -9 . 424291  47957 1 2E-0 i / ,  BR ( 6 )  / -3 . 5 1  12030408264E  +  00/ ,  BSE J  1  9  i  C 

BR  1 7 )  /-l. 572726362 0368E+ 01/,  BR ( 8  >  / -8. 228 1 439097 1 B6E+0 1 / ,  BSE J  1  9  1  D 

BR(9)  /-4. 9235537052367E+02/,  BR (10) /-3 . 3 1 621 85685480E  +  03/  BSE J 1 9 1 E 

BESJ 1 92 

DATA  C  1  v  1  >  /-2.0833333333333E-01/,  Cl (2)  /  1 . 2500000000000E-0 1  /  ,  BESJ193A 
Cl  (3)  /  0.0/,  Cl  (4)  /  0.0/,  Cl  (5)  /  0.0/,  Cl  (6)  /  0.0/,  BESJI93B 

Cl  (7)  /  0.0/,  C 1  ( 8 )  /  0.0/,  C 1  < 9 )  /  0.0/,  CK10)/  0.0/,  BESJ193C 

C  1  ( 1 1 )  /  0.0/,  Cl  ( 12) /  3. 34201 38888889E  —  oi / ,  BESJ193D 

Cl (13) /-4.0104166666667E-01/,  Cl(14)/  7 . 03 1 2500000000E-02/ ,  BESJ193E 
Cl  (15)/  0.0/,  Cl (16)/  0.0/,  Cl  ( 17) /  0.0/,  Cl(18)/  0.0/,  BESJ193F 

Cl  (19)/  0.0/,  Cl (20) /  0.0/,  Cl  (21)/  0.0/,  C1122)/  0.0/,  BESJ  1938 

Cl (23) /-1.0258125964506E  +  00/,  Cl(24)/  1 . 84646267361  1  1 E  +  0 0 / ,  BESJ193H 
Cl(25) / -8. 912109375 OOOOE -01/,  Cl(26>/  7 . 3242 1 87500000E-02/ ,  BESJ193I 
C 1  (27)  /  0.0/,  Cl  (28 ) /  0.0/,  Cl(29)/  0.0/,  Cl(30)/  0.0/,  BESJ193J 

Cl  (31)/  0.0/,  C 1 ( 32 ) /  0.0/,  C 1  ( 33 )  /  0.0/,  BESJ193K 

Cl  (34)/  4.6695844234262E  +  00/,  C 1  (35) /- 1 .  12070026 1 6223E+0 1 / ,  BESJ193L 
Cl  (36) /  8.7891 23535 1562 E+00/ ,  Cl (37) i-2. 3640869 1 4O625E  +  00/ ,  BESJ 1 9 3M 
Cl  (38)  /  1.  1215209960938E-01/,  Cl(39)/  0.0/,  Cl(40)/  0.0/,  BESJ193N 
Cl (41)/  0.0/,  Cl (42) /  0.0/,  Cl (43) /  0.0/,  Cl(44)/  0.0/,  BESJ1930 

Cl (45) /-2.8212072558200E+01/,  Cl(46)/  8. 46362 1 7674601 E+01 / ,  BESJ193P 
Cl (471/-9. 1818241543240E+01/,  Cl(48)/  4 . 253499B745388E+0 1 / ,  BESJ9318 
Cl (49)/-7.3687943594796£+00/,  Cl (50)/  2. 27108001 70898E-01 / ,  BESJ193R 
Cl  (51 )  /  0.0/,  Cl (52)  /  0.0/,  Cl  (53) /  0.0/,  Cl(54)/  0.0/  BESJ193S 

BESJ 194 

DATA  C 1  ( 55 ) /  0.0/,  Cl(56)/  2.  1257013003922E+02/,  BESJ195A 

Cl (57)/-7. 65252468 14118E+ 02/,  Cl(58)/  1 . 0599904525280E+03/ ,  BESJ195B 
Cl (59) /-6. 9957962737613E+02/ ,  Cl(60)/  2.  1819051  174421E  +  02/ ,  BESJ195C 
C1(61)/-2.6491430486952E+01/,  Cl(62)/  5. 7250142097473E-01 / ,  8ESJ195D 
Cl ( 63 ) /  0.0/,  Cl (64) /  0.0/,  Cl<65)/  0.0/,  Cl  ( 66 ) /  0.0/,  BESJ195E 

Cl  (67>/-1.9194576623184E*03/,  CK68)/  8. 061722181 7373E  +  03/ ,  BESJ195F 
Cl (69) /-l . 3586550006434E+04/,  Cl(70)/  1 . 1 655393336864E+04/ ,  BESJ195G 
Cl (71) / -5 . 3056469786 13 4E  +  03/ ,  Cl(72)/  1 . 2009029 1 32 1 64E  +  03 / ,  BESJ195H 
Cl <73>/-1.0809091978840E+02/,  CK74)/  1 . 7277275025845E+00/ ,  BESJ195I 
C 1  ( 75 )  /  0.0/,  Cl  ( 76 )  /  0.0/,  CH77)/  0.0/,  BESJ  1 95 J 

Cl ( 78 ) /  2. 02042913309 66E+04/,  Cl  (79) /-9. 69B0598388638E+04/ ,  BESJ 1 95K 
Cl (80)/  1 . 9254700123253E+05/ ,  Cl (81 ) / -2 . 03400 1 772804 2 E *05 / ,  BESJ195L 
Cl ( 82 ) /  1.2220046498302E+05/,  Cl (83) /-4.  1  1 92654968898E  +  04 / ,  BESJ195M 
Cl (84)/  7. 1 095 1 4302 489 4E +03/ ,  Cl ( 8 5 ) / - 4 . 9 3 9 1 5 304 7 7 309E +02 / ,  BESJ195N 

C 1  ( 86 )  /  6.0740420012735E+00/,  Cl(87)/  0.0/,  CK88)/  0.0/  BESJ1950 

BESJ 196 

DATA  C 2 ( 1 )  /-2. 42919187900 55E+05/,  C2(2)  /  1 . 3 1 1 7636 1 46630E+06/ ,  BESJ197A 
C 2 ( 3 )  / - 2 . 9980 1 59 1 8538 1 E  +  06 / ,  C2(4>  /  3.  76327  1  2976564E  +  06/ ,  BESJ197B 

C 2 ( 5 )  /-2.8135632265865E+06/,  C2(6)  /  1 . 2683652733216E  +  06/ ,  BESJ 1 9 7 C 

C 2 ( 7 )  /-3.3164517248456E  +  05/,  C2  (8)  /  4 . 52  1 876898  1  363E  +  04 / ,  BESJ 1 97D 

C2 ( 9 )  / -2, 49983048181  12E+03/,  C2(  10)/  2. 4380529699556E  +  0 1  /  ,  BESJ197E 

C 2 < 1 1  i  /  0.0/,  C2 ( 1 2) /  3. 284469853 0720E+06/,  BESJ197F 

C  2  < 13) /-l. 97068191  18432E  + 07/,  C  2 ( 1 4 ) /  5 . 0952602492665E  +  07  / ,  BESJ 1978 
C 2 ( 1 5 ) /-7. 41051 4821  1533E  +  07/ ,  C2 ( 1 6 ) /  6. 6344512274729E  +  07/ ,  BESJ  1 9 7 H 
C2(17)/-3. 7567176660763E+07/,  C2(18)/  1 . 3288767166422E+07/ ,  BESJ197I 


C2  ( 19) / - 2 . 7 8 5 6 1 8 i 2 8 0 8 6 4 E  +  0 6 / 
C2 ( 2 1 > / -  1 . 38860897537 17E+04/ 


DATA  A 1  ( 1 ) 

1  A  1  (  3 ) 

2  A  1  i  5  ) 


A  1  (  7  >  / 
A 1  ( 9 )  / 

A  1  ( 1  1 )  / 
A 1  <  1 3)  / 
A 1  (15)/ 
A 1  <  1 7 )  / 

A 1  (19)  / 
At  (21)  / 

A 1  (  23  )  / 
Al  (25)  / 
A1  (27)  / 
Al (29) /- 
A  1 (31) / - 
A 1  (33) / - 
Al  (35) /- 
Al  (37) /• 


4444444444444E-03/ 
848928B4892B8E-05/ 
4669 137274 179E -04/ 
61824297 06150E- 04/ 
3272104008323E-04/ 
0073B85876275E-04/ 
730607759 1788E- 04/ 
502747741609 1 E-04/ 
3166881654592E-04/ 
1640527147474E-04/ 
037724104 2299E -04/ 
3212051724950E-05/ 
429631057 1570E-05/ 
9373554 1 35459E-04/ 
4198627 3 55669E- 05/ 
5  0  8  0  3  5  5  8  0  5  3  0  5  E  -  0  4  / 
4680975664647E-04/ 
1974497568425E-04/ 
37699549891  19E-0 5/ 


DATA  Al (39) /• 
A  1  (  4  1  )  /  - 
l  A 1 ( 43 ) /  - 

S  A 1 ( 45 ) /  - 

)  A  1 (47) / - 

i  A  1  ( 49 )  /  ■ 

j  Al  (51)/- 


7. 2937434815522E-05/ 
5.6961156600937E-05/ 
4. 48 13486800888E-05/ 
3.5540053297204E-05/ 
2. 83996793904 18E-05/ 
2. 28459297 16472E- 05/ 
1.8481621762767E-05/ 


DATA  A2 ( 1 ) 


A2  (  3  )  / 
A2 (5 )  / 
A2  (  7  )  / 
A2 ( 9  )  / 
A2  ( 1 1 ) / 
A2( 13) / 
A2  ( 15) / 
A2 ( 1  7  )  / 
A2  ( 1 9 ) / 
A2  <  2 1 ) / - 
A2  ( 23 ) /- 
A 2  <  25  )  /  - 
A2  (27)  / 
A2(29)/- 
A2 (31 )  /- 
A  2 (33) / - 
A2 ( 35 )  /- 
A  2(37)/- 


.  5421  1 97 1A5774E— 04/ 
.  0446550359494E-05/ 
.  67471  1066997  IE— 04/ 
.  56501 4276 0860E -04/ 
.  1488669202982E-04/ 
,  644984 1925090E- 05/ 
,  743942992905 IE-05/ 
,  699397 1497922E-05/ 
.  3005667479396E-05/ 
,  59257 4858 1935E -06/ 
,  6539676969794E-06/ 
,  72739296 09125E- 06/ 
.  3130471509356E-06/ 
.  781941992 0177E- 04/ 
,  3793850631886E-05/ 
,  10916256 02736E- 04/ 
,  7895027379 132E- 04/ 
,  7524528 034085E-0 4/ 
,  2298287282 021E- 05/ 


DATA  A2  (39) /-1 . 7233430236696E-05/ 


A2 ( 4 1 ) /  2. 3139544314829E-05/ 
A 2  (43) /  4. 580061244 90 19E-05/ 
A2  ( 45 ) /  5. 6839620854582E-05/ 


C 2  (20) /  3. 0818640461 266E  +  05/ , 
C 2  (22) /  1.  1001714026925E+02/ 


A 1  ( 2 )  /- 
A  1  (  4  )  / 

A 1  ( 6  )  / 

A 1  ( 8  )  / 

A  1  (  1 0 )  / 

A 1  (12)/ 
Al  (14)/ 
Al ( 16) / 
Al ( 18) / 
Al (20) / 
At (22)  / 
Al (24) / 
Al (26)  / 
Al (28) / 
Al (30) 1- 
A 1 (32) /- 
Al (34) /- 
A 1  (36)/- 
A 1  (38)  /- 


.  2207772207792E- 
.  6592768783245E- 
.  6599558934626E- 
,  4873043734466E- 
.  1  63624857 1 236E- 
,  8626763663754E- 
,  6109170592902E- 
,  4050349739 1 27E- 
.  2366744559825E- 
.  0979829837271E- 
.  8262607836936E- 
,  857 1 08524787 1 E- 
.  0349754840779E- 
.  3224 1 745 1 8292E- 
,  1644493 1 67205E- 
,  551219249181 0E- 
,  3381550386749E- 
.  06 18431 920797E- 
,  26923045588 1 9E- 


A 1 (40)  / 
Al (42)  / 
Al (44)/ 
Al (46)/ 
Al (48)/ 
Al (50) / 
Al (52) / 


-6.4404235772102E-05/, 
-5.0473104430356E-05/, 
-3.9868872771760E-05/, 
-3. 1 74 1 425660902E— 05/ , 
-2.5452272063487E-05/, 
-2.0535275310648E-05/, 
-1. 665193300213 9E -05/ 


A2  ( 2 )  /- 
A2  ( 4 )  / 
A2 (6)  / 
A2  (8 )  / 
A2 ( 10) / 
A2 ( 1 2 )  / 
A  2  ( 1 4 )  / 
A2 ( 1 6 )  / 

A  2  <  1 8 )  / 

A  2  ( 2  0 )  / 
A 2 ( 2 2 )  /- 
A2 (24)  /  - 
A  2 (26) /- 
A2  ( 28 )  / 

A  2 (30) /- 
A2 (32) /- 
A2 ( 34 )  /- 
A2<36) /- 
A  2(38)/- 


.5616 12639451 6E- 
.  30 1 9865577324E- 
,  702225876B359E- 
,  36339 1 709774 4E- 
,  4586909303469E- 
.  0757033496520E- 
,  6275751 200534E- 
.  932 1093824794E- 
,  8262086674450E- 
440400498 1 425E- 
9134686709849E- 
,  1726937967866E- 
020114 1 879802E- 
,02471 95276 1 82E- 
,  3859823060301E- 
136801 l 524758E- 
,  2856408261 9 1 4 E - 
25544063 06069E- 
628607305881 2E- 


A2  ( 40 ) /  5. 6069048230460E-06/ , 


A2  (42) /  3. 6264274585679E-05/ , 
A2  ( 44 ) /  5. 245952949591  IE-05/ , 
A2  (46) /  5. 943498203931 0E-05/ , 
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BESJ199Q 
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BESJ199S 
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BESJ201A 

BESJ201B 

BESJ201C 

BESJ201D 

BESJ201E 
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BESJ202 

BESJ203A 

BESJ203B 

BESJ203C 

BESJ203D 

BESJ203E 

BESJ203F 

BESJ2036 

BESJ203H 

BESJ203I 

BESJ203J 
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BESJ2030 
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BESJ203Q 

BESJ203R 

BESJ203S 

BESJ204 

BESJ205A 

BESJ205B 

BESJ205C 

BESJ205D 


c. 


'*5' 


4  h  2  ( 4  7  i /  6.0647852757842E-05/ ,  A 2 ( 4 3 ) /  6. 0802390778844E-05/ , 

5  A 2  (49) /  6. 0157789453946E-05/ ,  A 2 (50) /  5 . 89 i 9965734470E-05/ , 

6  A  2  1 5  1 ! /  5.72515S2377759E -05/,  A2 ( 52  > /  5. 5280437558585E-05/ 

DATA  81  (  1)  /  1 . 7998872141355E-02/ .  81  (2)  /  5. 5996491  106439E-03/ , 

1  81  (3)  /  2. 88501402231  13E-03/ ,  81(4)  /  1 . 8009660676105E-03/ , 

2  81(5)  /  1. 247531 1 058920E-03/ .  81(6)  /  9 . 2287887657294E-04 / , 

3  81(7)  /  7. 1443042172729E-04/ ,  81(8)  /  5 . 7 1 78728 1 78970E-04 / , 

4  81(9)  /  4 . 6943100760648E-04/ ,  81(10)/  3. 9323283546292E-04/ , 

5  B 1  (11)/  3. 34818889 3 183 0E- 04/,  81(12)/  2 . 88952 1 4849575E-04 / , 

6  81(13)/  2. 5221  161554957E-04/,  81(14)/  2 . 2228058079888E-04  /  , 

7  B 1  ( 1 5 ) /  1. 9754 18380 33 06E- 04/,  81(16)/  1 . 76836B5501972E-04/ , 

8  81(17)/  1. 59- 3 168996618 2E -04/,  81  (  18)/  1 . 44347930 1 9733E-04  /  , 

9  81(19)/  1 . 314480681  1 996E-04  /  ,  81  (20)/  1 . 2024544494930E-04  /  , 

1  81  (21)/  1  .  1044914450460E-04/,  81(22)/  1 . 0 1 82877074057E-04  /  , 

2  B 1 (23) /  9, 419982242 0424E- 05  /  ,  Bl(24)/  8. 74  1 3054575383E-05/ , 

3  81(25)/  8. 1346626216280E-05/,  BK26)/  7 . 5900226964622E-05  / , 

4  B1 (27) /-l.  4928295 321343E-03/  ,  81  (28) /-B.  78204 70954639E-04 / , 

5  81  (29) / -5. 029 1654957 204E -04/  ,  81  (30) / -2 . 94822  1  385 1  275E-04 /  , 

6  81 (31  ) /-1  .  7546 3996970 78E- 04/  ,  81  ( 3 2 ) / -  1 . 0 4 0 0 8 5 5 0 4 6 0 8 2 E - 0 4 /  , 

7  81(33) /-5.9614195304646E-05/  ,  81  (34) / -3. 1 2038929076 1 0E-05/  , 

8  B1 (35) /-1 . 26089 7 35 98023 E-05/  ,  81  (36) / -2. 428926085 75 73E -07/  , 

9  81(37)/  8.0599616541427E-06/,  81(38)/  1 . 36507009262 l 5E-05/ 

DATA  81  (39)/  1 . 73964 1 2547293E-05/ ,  81  (40)/  1 . 98672978842 1 3E-05 / , 

1  B 1  (41  )/  2.  1446326379082E-05/,  BK42)/  2. 2395465923246E-05/ , 

2  81  (43)/  2.2896778381471E-05/,  81(44)/  2. 3078538981  1  1BE-05/ , 

3  81(45)/  2. 3032 1 9760809 1 E-05/  ,  Bl(46)/  2. 2823607372035E-05/ , 

4  81(47)/  2. 25005881  10529E-05/,  BK48)/  2 . 2098 10 1 536 1 99E-05/ , 

5  B 1 ( 49 ) /  2. 1641842744810 E-05/,  BK50)/  2 . 1 150764925622E-05/ , 

6  81(51)/  2.0638874978217E-05/,  81(52)/  2. 01 16524199708E-05/ 

DATA  82(1)  /  5.52213076721 29E-04/  ,  82  (2)  /  4 . 4793258  1  5523BE-04  /  , 

1  82(3)  /  2. 7952  0  653992  02E-04/',  82  (4)  /  1 . 52468 1 56 1 9845E-04  /  , 

2  82(5)  /  6. 9327 1  10565704E-05/' ,  82(6)  /  1 . 7625868306999E-05/ , 

3  82  (7)  /-1.3574499634327E-05/,  82  (8)  / -3 .  1 79724  1 335043E-05  /  , 

4  82(9)  /-4. 1 886 1 86 1 6 9669E-05/ ,  82 (  1 0) / -4 . 69004889379 1 4E-05/ , 

5  82 ( 1  1 > / -4 . 876654  4  741 379E-05 /  ,  82 ( 12) /-4. 870 1003 1  1 8674E-05/ , 

6  B2< 13) / — 4 . 74755620 890 09E-05/ ,  B2 ( 1 4 ) / -4 . 558 1 3058 1 3863E-05/ , 

7  82 ( 15) /-4. 3330964451  1 2 7 E - 0 5 /  ,  B 2  (  1 6 ) / - 4 .  0 9 2 3 0  1 9 3  1  5 7 7 5 E - 0 5 / , 

8  82(17) / -3. 348226 3 86 03 22E- 05/,  82 ( 1 8) / -3. 60857 1 675354 1 E-05/  , 

9  82  (  19) / - 3 . 3 7 7 9 3 3 0 6 1 2 3 3 7 E - 0 5 /  ,  82  (20) / -3.  1 5888560772 1  1 E-05  /  , 

1  82 (21 >/-2. 95269561 75081 E-05/ ,  82 ( 22 ) / -2 . 759789 1 482834E-05/  , 

2  82(23) / -2. 580061 7466683E-05/ ,  82 (24) /-2. 41 30835676 128E-05/ , 

3  82(25) / -2. 258235095 1 835E-05/  ,  82  (26) / -2.  1  1  4  796567689 1  E-05 / , 

4  B2  (27) / - 4 , 7461779655996E-04/  ,  B2  (28) / -4 .  77864567 1 4  732E-04 /  , 

5  82(29) / - 3 . 20390228 0 67 04E- 04/ ,  B2 (30) /-l . 61 1050161 1996E-04  /  , 

6  82 (31 i / - 4 . 257781 0128544E-05/ ,  82(32)/  3. 4457 1 29429497E-05/ , 

7  82  (33)/  7, 97092684 0  7568E- 05  /  ,  82  (34)/  1 . 03 1 3823670B27E-04  /  , 

9  82  (35)/  1  .  1246677526220E-04/ .  82  (36)/  l .  1 3 1 03642 1 0848E-04 / , 

9  82(37)/  1 . 0865163484877E-04/ ,  82(38)/  1 . 0  1  43795 1 59766E-04 / 

DATA  82  (39)/  9 . 2929839659336E-05/ ,  82(40)/  8 . 40293 1330 1 609E-05/ , 

1  82  (4  1)/  7. 52727991  3491  3  E-05/',  82(42)/  6. 6963252 1  97573E-05/ , 

2  82  (43)/  5. 9256454732320 E-05/,  82(44)/  5 . 22 1 6930882698E-05/ , 
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3  B  2  ( 4  5  > /  4. 58539485165 3 6E- 05/,  B2  (46)/  4 . 0 1 4455 1 389 1 49E-05 / , 

4  B2  ( 47 ) /  3. 5048 173003 133E-05/  ,  B2  ( 48 ) /  3. 0515799503435E-05/ , 

5  82  (49)/  2. 7  49561  1995052E-05/ ,  82(50)/  2. 2936363369 1 00E-05/ , 

6  82(51)/  1 . 9789305666402E-05/  ,  82(52) /  1 . 7009 1 9846364 1 E-05/ 

DATA  83(1)  /  7. 3646581 057258E- 04/,  83 ( 2 )  /  8. 7279080514619E-04/ , 

1  B3  ( 3 )  /  6. 2261486257314E  —  04  /  ,  83(4)  /  2 . 85998 1 54 1 9430E-04/ , 

2  B3(5)  /  3. 84737672879 3 7E- 06/,  B3(6)  / -  1 . 8790600363697E-04/ , 

3  83  (7)  /-2. 9760364659456E-04/  ,  83(8)  /-3. 4599812683267E-04/ , 

4  83(9)  /-3. 533824709 1604E-04/ ,  B3 ( 10) /-3. 3571563577505E-04/ , 

5  B3  (  1  1 ) /-3. 04321  124789 04E- 04  /  ,  B3  (  12)  1-2. 667227 230476 IE -04/ . 

6  B3(13)/-2. 276542141228 2E -04/,  B3 (  1 4 )/- 1 . 399226 1  1 85456E-04 / , 

7  B3< 15) /-1 . 550589 1859909E -04/,  B3  ( 1 6 >/- 1 . 23778240761 87E-04 /  , 

8  d3 (17)/ -9. 629261 4771 764E-05/,  83 ( 18) /-7. 251 7832771 442E-05/ , 

9  B3  (  19) /-5. 2207002889563E-05  /  ,  B3 (20) / -3. 5034775051 190E-05/ , 

1  B3  (21 ) /-2. 06489761 03555E- 05/,  B3 ( 22 ) / -8 . 70 1 0609684977E-06/ , 

2  B 3 (23)  /  1  .  1369868667510E-06/,  B 3  ( 2 4 ) /  9.  1642647412278E-06/ , 

3  63(25)/  1 . 5647778542887E-05/,  83(26)/  2. 082236294B247E-05/ 

DATA  GAMA (  1 )  /6 . 2996052494744E-0 1 / , 

1  GAMA  (  2 )  /2. 519842 0997898E- 01/,  GAMA (3)  / 1 . 54  7903004 1 566E-0 1 / , 

2  GAMA ( 4 )  /l. 1071306241616E-01/,  GAMA ( 5 )  /8. 5730939552740E-02/ , 

3  GAMA ( 6 )  /6.9716131695868E-02/,  GAMA  ( 7 )  /5. 8608567 1 3937 1 E-02/ , 

4  GAMA  ( 8 )  /5.0469887353631E-02/,  3AMA(9)  / 4 . 42600580689 1 6E-02/ , 

5  GAMA( 10) /3. 9372066 154351 E-02  / ,  GAMA (  1  1) /3.  54283 1 9 59 2 4 4 6 E - 0 2 / , 

6  GAMA ( 12) /3. 2 18 18857502 10E-02/,  GAMA ( 1 3 ) / 2 . 9464624079 1  1 6E-02/  , 

7  GAMA ( 14) /2. 715816771 1293E-02/,  GAMA (15) /2.5176S272973B6E-02/, 

8  GAMA ( 16) /2. 345707 55 3060SE -02/ ,  GAMA (  17) /2.  1 9508390 1 349 1 E-02/  , 

9  GAMA (18) /2.0621082823565E-02/,  GAMA ( 19) / 1 . 9438824089738E-02/  , 

1  GAMA (20) /  1.838 1 063380068E-02/ ,  GAMA ( 2 1 ) /  1 . 742932 1 323196E-02/ , 

2  GAMA (22)/  1.6568583778661  E-02/,  GAMA (23 ) / 1 . 573652B59B792E-02/ , 

3  GAMA (24) / 1 . 5072950 1494 10E-02/ ,  GAMA (25) / 1 . 4419325083996E-02/ , 

4  GAMA ( 26 ) / 1 . 38 1 8480573534E -02 / 


TEST  INPUT  ARGUMENTS 

NZ=0 
KT=  1 

IF(N-l)  92,108,109 

108  KT  =  2 

109  NN=N 

I F  <  X )  93,110,120 

110  I F (ALPHA )  91,114,116 
114  Y ( 1 ) = 1 . 

IF(N.EQ.i)  RETURN 
11=2 

GO  TO  118 
116  11=1 

118  DO  119  I  =  1 1  ,  N 

119  Y(I)=0. 

RETURN 

120  CONTINUE 

IF (ALPHA. LT. 0. )  GO  TO  91 
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DFN=DBLE (FLOAT ( N) >+DBLE( ALPHA) -1  .  D+O 

FNU=DFN 

X02=X*. 5 

SX02=X02*X02 


DECISION  TREE  FOR  REGION  WHERE  SERIES,  ASYMPTOTIC  EXPANSION  FOR 
TO  INFINITY  AND  ASYMPTOTIC  EXPANSION  FOR  NU  TO  INFINITY  ARE 
APPLIED. 


IF  (SXQ2. LE . (FNU  + 1 . ) )  GO  TO  850 
TA=AMAX1 (20. ,FNU> 

IF(X.GT.TA)  GO  TO  BBO 
IF  ( X . GT . 12. )  GO  TO  860 
X02L=AL0G ( XQ2) 

NS=SX02-FNU 
60  TO  852 
850  FN=FNU 

FNP 1 =FN+ 1 . 

X02L=AL06<X02) 

I  S  =  KT 

IF (X . LE. 0. 5)  GO  TO  134 
NS=0 

852  DFN=DFN+DBLE (FLOAT (NS) ) 

FN=DFN 

FNP1=FN+1. 

IS=KT 

IF (N-l+NS. GT. 0)  I S=3 
GO  TO  134 

860  NS=AMAX1 <36.-FNU,0. ) 
DFN=DFN+DBLE (FLOAT (NS) ) 

FN=DFN 


I S  =  K  T 

IF  (N-l+NS. GT. 0)  I S=3 


GO  TO  130 
880  CONTINUE 
RTX  =  SQRT  ( X ) 

TAU=RTWO*RTX 
TA=TAU+FNULIM (KT) 
IF(FNU.LE.TA)  GO  TO  500 
129  FN=FNU 
I  S=KT 


UNIFORM  ASYMPTOTIC  EXPANSION  FOR  NU  TO  INFINITY 


130  CONTINUE 
X  X  =  X  /  F  N 
W2= 1 . -X  X *X  X 
ABN2  =  ABS ( W  2  > 

RA=SQRT ( ABW2 ) 

IF(ABW2.GT. 0.2775)  GO  TO  200 


CASES  NEAR  X  =  FN,  ABS ( 1 . - ( X /FN ) **2 ) . LE . 0 . 2775 
COEFFICIENTS  OF  ASYMPTOTIC  EXPANSION  BY  SERIES 
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ZETA  AND  TRUNCATION  FOR  A(ZETA)  AND  B(ZETA)  SERIES 

KMAX  IS  TRUNCATION  INDEX  FOR  A (ZETA)  AND  B(ZETA)  SER I ES  =  MAX  ( 2 , SA ) 

S  A  =  0 . 

IF ( ABW2. EQ. 0. )  GO  TO  21 
SA=TOLS/ALOG (ABW2) 

21  S8=SA 

DO  22  1=1,5 

KMAX (I ) =AMAX 1 (SA,2. ) 

SA=SA+SB 

22  CONTINUE 
KB=KMAX (5) 

KLAST=KB- 1 
SA=GAMA (KB ) 

DO  24  K= 1 , KLAST 
KB=KB- 1 

SA=SA*W2+GAMA(KB) 

24  CONTINUE 
Z=N2*SA 
A  Z  =  A  B  S  (  Z  ) 

RTZ=SQRT ( AZ ) 

FN13=FN»*C0N2 

RTARY=RTZ*FN1 3 

ARY=-RTARY*RTARY 

AZ32=AZ*RTZ*C0N 1 

ACZ=FN*AZ32 

IF  (Z . LE. 0. )  GO  TO  27 

TEST  FOR  UNDERFLOW,  l . E-280=EXP (-644. ) ,  ONE  WORD  LENGTH 
UP  FROM  UNDERFLOW  LIMIT  OF  CDC  6600 

IF  (ACZ. GT. ELIM2)  GO  TO  180 
ARY=-ARY 

27  PHI=SQRT (SORT (SA+SA+SA+SA) ) 

B  ( ZETA)  FOR  S  =  0 

KB=KMAX (5) 

KLAST = K.  B  —  1 
SB=BETA(KB, 1) 

DO  23  K=l, KLAST 
KB=KB- 1 

SB=SB*W2  +  BETA(KB,  1) 

23  CONTINUE 
K  S  P  l  =  1 
FN2=FN*FN 
RFN2= 1 . /FN2 
RDEN= 1 . 

ASUM=  l . 

RELB=TOL*ABS (SB) 

8SUM=SB 
DO  25  KS=  1 ,4 
KSP 1 =KSP 1 + 1 
RDEN=RDEN*RFN2 
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C  A(ZETA)  AND  B(ZETA)  FOR  S= 1 .2,3,4 
C 

KB=KMAX (5-KS) 

KLAST=KB-1 

SA=ALFA(KB,KS) 

SB=BETA(KB,KSP1> 

DO  26  K= 1 , KLAST 
KB=KB-1 

SA=SA*W2+ALFA(KB,K5i 

SB=SB*W2*BETA(K8,KSP1) 

26  CONTINUE 
TA=SA*RDEN 
TB=SB*RDEN 
ASUM=ASUM+T  A 
BSUM=BSUM+TB 

I F ( ABS  <  TA ) .LE.TOL. AND. ABS(TB) .LE. RELB)  GO  TO  152 
25  CONTINUE 
152  CONTINUE 

BSUM=BSUM/ (FN*FN13) 

GO  TO  400 
C 

200  CONTINUE 
T  AU= 1 . /RA 
T  2  =  1 . /W2 

IF (W2. BE. 0. )  GO  TO  30 
C 

C  CASES  FOR  (X/FN) .GT.SQRT (1.2775) 

C 

AZ32=A8S(RA-ATAN(RA1) 

ACZ=AZ32*FN 
CZ=-ACZ 
Z32=l . 5*AZ32 
RTZ=Z32**C0N2 
FN13=FN**C0N2 
RTARY  =  RTZ*FN1  3 
ARY=-RTARY*RTARY 
GO  TO  150 
30  CONTINUE 
C 

C  CASES  FOR  (X/FN). LT.SQRT(0. 7225) 

C 

AZ32=ABS(AL0G< (l.+RA) /XX)  -RA) 

C 

C  TEST  FOR  UNDERFLOW,  l.E-280  =  EXP<-644,),  ONE  WORD  LENGTH 
C  UP  FROM  UNDERFLOW  LIMIT  OF  CDC  6600 
C 

ACZ=AZ32*FN 


BESJ350 

BESJ351 

BESJ352 

BESJ353 

BESJ354 

BESJ355 

BESJ356 

BESJ357 

BESJ358 

BESJ359 

BESJ360 

BESJ361 

BESJ362 

BESJ363 

BESJ364 

BESJ365 

BESJ366 

BESJ367 

BESJ368 

BESJ369 

BESJ370 

BESJ371 

BESJ372 

BESJ373 

BESJ374 

BESJ375 

BESJ376 

BESJ377 

BESJ378 

BESJ379 

BESJ380 

BESJ381 

BESJ382 

BESJ383 

BESJ384 

BESJ385 

BESJ386 

BESJ387 

BESJ388 

BESJ389 

BESJ390 

BESJ391 

BESJ392 

BESJ393 

BESJ394 

BESJ395 

BESJ396 

BESJ397 


150  CONTINUE 

PHI=SQRT( ( RTZ+RTZ ) *TAU> 

TB=  1 . 

ASUM=1. 

TFN=TAU/FN 

UPOL  (2)  =  (C ( 1 , 1 ) *T2  +  C (2, 1 ) ) *TFN 

RCZ=C0N1 / C  Z 

CRZ32=C0N548*RCZ 

BSUM=UP0L(2)+CRZ32 

REL8=T0L*ABS(BSUM) 

AP=TFN 

KS=0 

KP1  =  2 

RZDEN=RCZ 

DO  155  LR=2 , 9 , 2 


COMPUTE  TWO  U  POLYNOMIALS  FOR  NEXT  A(ZETA)  AND  B(ZETA) 


LRP 1 =LR+ 1 
DO  101  K=LR , LRP 1 
KS=KS+ 1 
KP1=KP1+1 
S 1 =C ( 1 , K ) 

DO  102  J=2,KP1 
S 1=S 1 *T2+C ( J , K ) 

102  CONTINUE 
AP=AP*TFN 
UPOL ( KP 1 ) =AP*S 1 
CR(KS)=BR(KS)#RZDEN 
RZDEN=RZDEN*RCZ 
DR(KS)=AR(KS)*RZDEN 
101  CONTINUE 

SUMA=UPOL ( LRP 1 ) 


SUMB=UP0L(LR+2)+UP0L(LRPl)*CRZ: 


J  U=LRP 1 

DO  151  JR=1,LR 
JU= JU- 1 

SUMA=SUMA+CR(JR)*UPOL(JU) 

SUMB=SUMB+DR(JR)*UPDL(JU) 

151  CONTINUE 
TB=-TB 

IF  (W2.GT.0. )  TB=ABS  (TB) 

ASUM=ASUM+SUMA*TB 

BSUM=BSUM*SUM8*TB 

I F  ( ABS ( SUMA ) .LE.TOL. AND. ABS(SUMB) .LE.RELB) 


SO  TO  165 


155  CONTINUE 
165  TB=RTARY 

IF (W2.GT. 0. ) 
BSUM=BSUM/TB 


400  CONTINUE 

CALL  JAIRY(ARY,RTARY,ACZ,AI,DAI) 
TEMPI  IS) =PHI*(AI*ASUM+DAI*BSUM)/FN13 
60  TO  (401,202,650)  ,  IS 
402  TEMP ( 1 ) =TEMP (3) 
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KT  =  1 

BESJ460 

401 

I  S  =  2 

BESJ461 

DFN=DFN-1 . D+0 

BESJ462 

FN=DFN 

BESJ463 

GO  TO  130 

BESJ464 

BESJ465 

SERIES  FOR  (X/2) **2. LE. NU+1 

BESJ466 

BESJ467 

134 

CONTINUE 

BESJ468 

GLN=GAMLN  (FNP 1 ) 

BESJ469 

ARG=FN*X02L-GLN 

BESJ470 

IF1ARG.LT. -ELIM1)  GO  TO  123 

BESJ471 

EARG=EXP1ARG) 

BESJ472 

300 

CONTINUE 

BESJ473 

S  =  l. 

BESJ474 

AK  =  3 . 

BESJ475 

T2  =  l . 

BESJ476 

T=  1 . 

BESJ477 

S 1  =FN 

BESJ478 

DO  125  K  =  1 , 1  7 

BESJ479 

S2=T2+S 1 

BESJ480 

T  =  -T  *SXQ2/S2 

BESJ481 

S=S  +  T 

BESJ482 

IF1ABS1T) .LT.TOL)  GO  TO  127 

BESJ483 

T2=T2+AK 

BESJ4B4 

A  K  =  A  K  +  2 . 

BESJ485 

S1=S1+FN 

BESJ486 

125 

CONTINUE 

BESJ487 

127 

CONTINUE 

BESJ4BB 

TEMPI  IS)=5*EARG 

BESJ439 

GO  TO  (301,202,600) ,  IS 

BESJ490 

301 

EARG=EAR6*FN/X02 

BESJ491 

DFN=DFN-1 . D+0 

BESJ492 

FN=DFN 

BESJ493 

IS=2 

BESJ494 

GO  TO  300 

BESJ495 

BESJ496 

SET  UNDERFLOW  VALUE  AND  UPDATE  PARAMETERS 

BESJ497 

BESJ498 

180 

Y (NN) =0. 

BESJ499 

NN=NN- 1 

BESJ500 

DFN=DFN-1 . D+0 

BESJ501 

FN=DFN 

BESJ502 

IF  ( N  N  —  1 )  170,171,130 

BESJ  503 

171 

KT=2 

BESJ504 

I S  =  2 

BESJ505 

GO  TO  130 

BESJ506 

123 

Y (NN) =0. 

BESJ507 

NN=NN- 1 

BESJ508 

FNP 1 =FN 

BESJ509 

DFN=DFN-1 . D+0 

BESJ510 

FN=DFN 

BESJ51 1 

IF(NN-l)  170,172,173 

BESJ512 

172 

KT  =  2 

BESJ513 

IS=2 

BESJ514 
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IF  (SX02. LE. FNP1 )  GO  TO  133 
GO  TO  130 

ARG=ARG-X02L+AL0G(FNP1> 

I F  ( ARG. LT . -EL  I M 1 )  GO  TO  123 

60  TO  134 

NZ=N-NN 

RETURN 


BACKWARD  RECURSION  SECTION 


CONTINUE 

NZ=N-NN 

IFtKT.EQ.2)  GO  TO  250 
CONTINUE 

BACKWARD  RECUR  FROM  INDEX  ALPHA+NN-1  TO  ALPHA 
Y (NN ) =TEMP  ( 1 ) 

Y (NN-1 ) =  TEMP (2) 

IF (NN. EQ. 2)  RETURN 
DX  =  X 

TRX=2. D+O/DX 
DTM=DFN*TRX 
TM=DTM 
K=NN+1 

DO  230  1=3, NN 
K  =  K-  1 

Y  (K  —  2) =TM*Y  <K— 1 ) — Y (K) 

DTM=DTM-TRX 

TM=DTM 

CONTINUE 

RETURN 

Y ( 1 ) =TEMP  ( 2 ) 

RETURN 


ASYMPTOTIC  EXPANSION  FOR  X  TO  INFINITY  WITH  FORWARD  RECURSION  IN 
OSCILLATORY  REGION  X.GT.MAX(20,  NU) ,  PROVIDED  THE  LAST  MEMBER 
OF  THE  SEQUENCE  IS  ALSO  IN  THE  REGION. 


CONTINUE 
I N=ALPH A-T  AU  +  2 . 

IF ( IN. LE. 0)  GO  TO  502 
1NP1  =  1N+  1 

DALPHA=ALPHA-FLOAT  ( INP1 ) 
KT  =  1 

GO  TO  511 
DALPHA=ALPHA 
IN=0 
I  S  =  KT 

ARG=X-PIDT*DALPHA-PDF 
SA=S I N ( ARB ) 

SB=COS ( ARG ) 

RA=RTTP/RTX 
£TX=8. *X 
DX=DALPHA 
DX  =  DX  +  DX 
DTM=DX*DX 
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I 

I  T2  =  DTM- 1 . D+0 

T2=T2/ETX 
S2=T2 

RELB=TOL*ABS (T2) 

T 1  =ETX 
Sl  =  l. 

FN=  1 . 

AK  =  8 . 

DO  504  K  = 1 , 1 3 
T1  =  T  1  +  ETX 
FN=FN+-AK 

(  DX  =  FN 

|  TRX=DTM-DX 

AP=TRX 

T2=-T2*AP/T  1 
S1=S1+T2 
T 1=T 1 +ETX 
AK=AK+8. 

FN=FN+AK 

-  D  X=FN 

!  TRX=DTM-DX 

AP=TRX 

1  T2=  T2*AP/T 1 

,  S2=S2+T2 

f  IF (ABS (T2) .LE.RELB)  GO  TO  505 

\  AK=AK+8. 

’  504  CONTINUE 

(  505  TEMP(IS)=RA*(S1*SB-S2*SA) 

GO  TO  (506,507)  ,  IS 

506  DALPHA=DALPHA+1. 

IS=2 

TB=SA 

5A=-SB 

i  SB=TB 

GO  TO  503 

FORWARD  RECURSION  SECTION 

507  IF (KT. EQ. 2)  GO  TO  250 

51  =  TEMP ( 1 ) 

S2=TEMP ( 2 ) 

T  X  =  2 . /X 
TM=DALPHA*TX 
IF(IN.EQ.O)  GO  TO  520 

FORWARD  RECUR  TO  INDEX  ALPHA 

DO  510  1  =  1,  IN 
S  =  S2 

52  =  TM*S2-S  1 
TM=TM+TX 
S1=S 

510  CONTINUE 

IF  (NN. EQ. 1 )  GO  TO  535 
S  =  S2 


BESJ570 

BESJ571 

BESJ572 
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BESJ574 

BESJ575 

BESJ576 

BESJ577 
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BESJ582 
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BESJ600 

BESJ601 

BESJ602 

BESJ603 

BESJ604 

BESJ605 

BESJ606 

BESJ607 

BE3J608 

BESJ609 

BESJ610 

BESJ61 1 

BESJ612 

BESJ613 

BESJ614 

BESJ615 

BESJ616 

BESJ617 

BE5J618 

BESJ619 

BESJ620 

BESJ621 

BESJ622 

BESJ623 

BESJ624 
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S2=TM*S2-S1 

BESJ625 

TM  =  TM-*-TX 

BESJ626 

S1=S 

PESJ627 

520 

CONTINUE 

BESJ62B 

C 

BESJ629 

C 

FORWARD  RECUR  FROM  INDEX  ALPHA  TO  ALPHA+N- 1 

BESJ630 

c 

BESJ631 

Y  ( 1 ) =S1 

BESJ632 

Y  (2) =S2 

BESJ633 

IF  (NN. EQ. 2)  RETURN 

BESJ634 

DO  530  1=3, NN 

BESJ635 

Y  ( I  >  =  TM*Y ( I  —  1 ) — Y ( 1—2) 

BESJ636 

TM=TM+TX 

BESJ637 

530 

CONTINUE 

BESJ638 

RETURN 

BESJ639 

535 

Y  ( 1 1 =S2 

BESJ640 

RETURN 

BESJ641 

C 

BESJ642 

C 

BACKWARD  RECURSION  WITH  NORMALIZATION  BY 

BESJ643 

C 

ASYMPTOTIC  EXPANSION  FOR  NU  TO  INFINITY  OR  POWER  SERIES. 

BESJ644 

C 

BESJ645 

600 

CONTINUE 

BESJ646 

C 

COMPUTATION  OF  LAST  ORDER  FOR  SERIES  NORMALIZATION 

BESJ647 

KM=AMAX 1(3. -FN , 0 . ) 

BESJ648 

TFN=FN+FLOAT (KM) 

BESJ649 

TA=(SLN  +  TFN-0. 9189385332-0. 0833333333/TFN) / (TFN+O.  5) 

BESJ650 

TA=X02L-TA 

BESJ651 

TB=-<1.-1.5/TFN)/TFN 

BESJ652 

IN=CE/ (-TA  +  SQRT (TA*TA-CE*TB)  )  +  l .  5 

BESJ653 

IN= IN+KM 

BES-J654 

60  TO  603 

BESJ655 

650 

CONTINUE 

BESJ656 

C 

COMPUTATION  OF  LAST  ORDER  FOR  ASYMPTOTIC  EXPANSION  NORMALIZATION 

BESJ657 

GLN=AZ32+RA 

BESJ658 

IF  (ARY. GT. 30. )  GO  TO  675 

BESJ659 

RDEN=(PP(4)*ARY+PP(3i ) *ARY+ 1 . 

BE5J660 

RZDEN=PP ( 1 ) +  PP  (2) *ARY 

BESJ661 

TA=RZDEN/RDEN 

BESJ662 

IF  (W2. LT.0. 10)  GO  TO  651 

BES J 663 

TB=GLN/RTARY 

BESJ664 

GO  TO  677 

BESJ665 

651 

TB=< 1.25992 1049+0. 1679B94730*W2) /FN13 

BESJ666 

GO  TO  677 

BESJ667 

675 

CONTINUE 

BESJ668 

TA=CQN1*TCE/ACZ 

BESJ669 

TA= ( (0. 0493827 160*T A -0. llllllllll) *TA+0. 6666666667) *TA*ARY 

BESJ670 

IFIW2.LT. 0. 10)  GO  TO  651 

BESJ671 

T6=GLN/RT  ARY 

BESJ672 

677 

IN=TA/TB+1 . 5 

BESJ673 

IF! IN.GT. INLIM)  GO  TO  402 

BESJ674 

603 

DX=FLOAT (IN) 

BESJ675 

DTM=DFN+DX 

BESJ676 

D  X  =  X 

BESJ677 

TRX=2 . D+O/DX 

BESJ678 

DTM=DTM*TRX 

BESJ679 
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TM=DTM 

TA=0. 

TB=TOL 
KK=  1 

605  CONTINUE 
C 

C  BACKWARD  RECUR  UNINDEXED 

C 

DO  601  1  =  1,  IN 
S  =  TB 

TB=TM*TB-TA 
TA  =  S 

DTM=DTM-TRX 

TM=DTM 

601  CONTINUE 

C  NORMALIZATION 

IF  ( KK . NE .  1)  GO  TO  604 
TA= (TA/TB) *TEMP (3) 
TB=TEMP (3) 

KK  =  2 
IN=NS 

IF(NS.NE.O)  GO  TO  605 
604  Y (NN) =TB 
615  NZ=N-NN 

IF  CNN. EQ. 1 )  RETURN 
S  =  TB 

TB=TM*TB-Tfi 

TA=S 

DTM=DTM-TRX 
TM=DTM 
K=NN- 1 
Y (K ) =TB 

IF (NN. EQ. 2)  RETURN 
KM=K- 1 
C 

C  BACKWARD  RECUR  INDEXED 
C 

DO  602  1=1, KM 
Y  (K-l ) =TM*Y (K) —  Y (K  +  l ) 

DTM  =  DTM-TR  X 
TM=DTM 
K  =  K  - 1 

602  CONTINUE 
RETURN 

C 

C 

C 

91  CONTINUE 
NZ=-2 
RETURN 

92  CONTINUE 
N  Z  =  -  3 
RETURN 

93  CONTINUE 
NZ  =  - 1 
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RETURN 

END 

SUBROUTINE  J A  I  RY  (  X , RX , C , A I , DA  I ) 


CDC  6600  ROUTINE  1-2-74 

JAIRY  COMPUTES  THE  AIRY  FUNCTION  A I ( X ) 

AND  ITS  DERIVATIVE  DAI  ( X )  FOR  JBESS 

INPUT:  X  -  ARGUMENT,  COMPUTED  BY  JBESS,  X  UNRESTRICTED 

RX  -  RX=SQRT ( ABS ( X ) ) ,  COMPUTED  BY  JBESS 
C  -  C=2.*(ABS(X>**1.5>/3. ,  COMPUTED  BY  JBESS 
OUTPUT:  A I  -  VALUE  OF  FUNCTION  AI(X) 

DAI  -  VALUE  OF  THE  DERIVATIVE  DAI ( X ) 

WRITTEN  BY  D.E,  AMOS,  S.L.  DANIEL  it  M.K  WESTON 


DIMENSION  AK1 (14) , AK2 (23) ,AK3(14) 

DIMENSION  AJP (1 9) , AJN ( 1 9) ,  A  1 1 5) ,B(15) 
DIMENSION  DAK  1  (14) , DAK2 ( 24 ) ,DAK3(14) 
DIMENSION  DAJP(19) ,DAJN(19) ,DA!15) ,DB(15) 
DATA  N1  ,N2,N3,N4/14,23, 19, 15/ 

DATA  Ml, M2, M3, M4/ 12, 21, 17, 13/ 

DATA  FPI 12, CONI ,C0N2,C0N3,C0N4 ,C0N5/ 

1  1 . 3089969389958E+00,  6 . 6666666666667E-0  1 , 

2  3. 8000458986729E-01 ,  8, 3333333333333E-01 , 
DATA  AKl(t)  /  2.2042309098779E-01/, 

1  AK1(2)  /-1.2529024278770E-01/,  AK1 (3)  /  1. 

2  AK 1(4)  /  8.2284415200634E-04/,  AKl (5)  /-2. 

3  AKl  ( 6 )  /  1. 6382428 017212E-05/,  AK  1  (  7 )  /  3. 

4  AKl (8)  /- 1 . 29621 99935933E-07/ ,  AKl (9)  /  8. 

5  AKl  (10)/  1. 539639686233 OE-11/,  AKK1D/-3. 

6  AKl  (12)/  2.0325325742363E-12/,  A K 1 < 1 3 ) / -  1 . 

7  AKl ( 14) /-5. 1616949778508E— 15/ 

DATA  A K 2 ( 1 )  /  2 . 74366 1 5086960E-0 1 / , 

1  A  K  2  ( 2 )  /  5. 397909 6973690E-03/,  AK2(3)  /- 1. 

2  A K 2 ( 4 )  /  4. 274275282487 5E -04/,  AK2(5)  / - 1 . 

3  AK2  ( 6 )  /  2. 88763 1 7 1 3 1 890E-05/ ,  AK2(7)  1-1. 

4  AK2  ( 8 )  /  l .  8729020.974 102E-06  /  ,  AK2(9>  /-4. 

5  AK2 ( 1 0 ) /  1 . 21 13041695591E— 07/ ,  AK2(ll>/-3. 

6  AK2  (12)/  7.9245470523265E-09/,  AK2  ( 1 3 > / -2 . 

7  A  K  2  (  1 4 ) /  5. 2686305659574E-10/ ,  AK2(15)/-1. 

9  AK2 (16)/  3. 5614103 90 13 7! E- 11/,  AK2(17)/-9. 

9  AK2  (  18) /  2. 4446445 047364E-12/,  AK2(19)/-6. 

1  AK2  (20) /  1 . 70106 03055935E-13/ ,  AK2  <  2 1 ) / -4 . 

2  A  K  2  (22) /  1.  1977479916481E-14/,  AK2(23)/-3. 
DATA  AK3 ( 1 )  /  2 . 8027 1 44734079E-0 1 / , 

1  AK3  ( 2 )  /~1  .  781270 4284438E-03/,  AK3(3)  /  4. 

2  A K 3 ( 4 )  /-1.6324996526900E-06/,  AK3(5)  /  9. 

3  A K 3  ( 6 )  /-6. 5229433 022916E- 09/,  AK3(7)  /  5, 

4  AK3  (8)  /-5. 244082518 0026E-11/,  AK3(9)  /  5. 

5  AK3(  10) /  —  6. 5637524463931E— 13/ ,  AK3  ( 1 1 ) /  8. 

6  AK3  (  12) / - 1 . 1270513469106E— 14/ ,  AK3  (13)/  1. 

7  AK3< 14) /-2. 4648032431 243E-16/ 

DATA  AJP ( 1 !  /  7. 789529664 3 75SE-02/, 

1  AJP ( 2 )  /-1 . 8435636345680E-01/ ,  AJP  ( 3 )  /  3. 


5.0315471619678E+00, 

8.6602540378444E-01/ 

0388 1 1 63359 1 9E-02/ , 
3461 434589 123E-04  /  , 
06902589573 1 9E-07 / , 
22908 15882367E-09/, 
3916546561568E-11/, 
1067954609788E-14/, 


57339220621 19E-03/, 
1212491 739992E— 04  /  , 
3680422537055E-06/  , 
7589279396229E-07/  , 
092453742706 1 E-08/  , 
03902447 1 679 IE- 09  /  , 
3670476763957E-10/  , 
3138829654843E-12/, 
4384026199 096E-13/, 
507601 04 50328E- 14/, 
19077040865 07E-15/ 

0342257962900E-05  /  , 
2118143247677E-08/, 
471384 0457655E-10/  , 
604779041  1 72 1 E- 1 2/  . 
3128576196625E— 14/  , 
622679765981 3 E-15/, 


0141260521617E-02/, 
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3  AJN (6)  /  -  5 

4  A  J  N  ( 8  )  /  2 


DATA  BID 


2  A J F  >  4 )  /  3 . 053427242776  IE-02/ ,  AJP(5)  /-4. 

3  AJPto)  /-i . 7274955256395E-03/ ,  AJP(7>  /  2. 

4  A J  P  (  8 )  /  5. 0456477751708E-05/ ,  AJP(9)  /- 6. 

5  A  J  P  ( 10)  /  -  9 . 0398674551 077E  - 07  /  ,  AJP<11)/  9. 

6  A  J  P  ( 1 2 ) i  1 . 0963945330520E-08/ ,  AJP(13)/-l. 

7  AJP  ( 1  4)  /  -  9 . 6  035944 1 34465E- 1 17  f  A  J  P  ( 1 5  > /  8. 

8  AJP  (  16) /  6.3612343901877E-13/,  AJP(17)/-4. 

9  AJP  <  19)  /-3.2981028892962E-15/,  AJ P  (  1 9  )  /  2. 
DATA  A  J  N  (  1  )  /  3.8049788761724E-02/ , 

1  A  J  N  (  2  )  /-2.4531954194555E-01/,  AJN<3;  /  1. 

2  A  J  N  <  4 )  /  7. 4  93300458 1879E-02/ ,  A JN ( 5  >  /- 2. 

3  AJN  (6)  / -  5. 92535597 3 0498E- 03/,  AJN  (  7 )  /  1. 

4  AJN  (8)  /  2,  1831  183132222E-04/,  A JN ( 9 >  /-4. 

5  AJN  1 10)  /  -4. 66874994  1 7 177E  —  06/ ,  AJN(ll)/  7. 

6  AJN (  12) /  6. 529647708546 3E- 08/,  AJN  ( 13) / - 8 . 

7  AJN( 14) /-to. 441S615897698E-10/ ,  A  J  N  < 1 5 ) /  7. 

8  A  J  N ( 1 6 ) /  4. 724654 3 171785E-12/,  AJN(17)/-4. 

9  AJN  ( 18) / -  2 . 6776271 038919E-14/ ,  AJN119)/  2. 

DATA  Ail)  /  4.9027542474279E-01/ ,  A (2)  /  1. 

1  A  (  3 )  /  - 9 . 6619596 3 14 03  IE- 05/ ,  A(4)  /  1. 

2  A  ( 5 )  /  2.9815734265486E-07/,  A ( 6 )  /-l. 

3  A  ( 7 )  /-1.0368573766714E-09/,  A(8)  /  3. 

4  A  (  9 )  /-2. 5709141063278E-1  1/  ,  AU01/-2. 

5  A  ( 1  1  ) /  9. 57523279 04826E-1 3/,  A  ( 1  2 ) / - 1 . 

6  A  ( 1 3 ) / -2. 9090771677 072E- 15/ ,  A  i  1  4 ) /  4. 

7  A ( 15) / -9 . 990038748 1 026E- 1 6 / 

DATA  Bil)  I  2. 7859355280 30 8E-01/ ,  8(2)  /~3. 

1  B  (  3 )  /-2.31  14967738499E-05/,  B ( 4 )  /  4. 

2  6(5)  /-l . 1 241 590793 1 33E-07/ ,  B(6)  /-2. 

3  8(7)  /  2. 6 0948 0753 0219c- 09/,  8(B)  /-3. 

4  B (9 )  /-3.5084997842388E-11/,  8(10)/  5. 

5  B(lli/-2.0464482875333E-13/,  B  ( 12) /-I . 

6  8(13)/  2. 37724778 0 3 878E-14/,  B(14)/-2. 

7  B(15> /- 3. 3265631 169617E-16/ 

DATA  N1D,N2D,N3D.N4D/14,24, 19, 15/ 

DATA  M1D,N2D,M3D,M4D/12,22, 17,  13/ 

DATA  DAK  1(1)  /  2. 0456784230789E-0 1 / , 

1  DAKi(2)  /  -  a . 6 1 3227 399056 6E -02 / ,  DAK1(3)  /- 

2  DAKK4)  /  3.  1  2  1 8 3 4 9  1  5 5 6 2 9 E - 03 / ,  DAK1(5)  /- 

3  DAK  1  (6)  /-6.  35636298679 3 9E- 06/ ,  D A K  1  ( 7 )  / 

4  DAK  1(8)  / -2 .  1831  1  1 9533009E-07/ ,  D A K  1  ( 9 )  /- 

5  DAK  1(10)/  1 .  130980356223  IE -09/ ,  DAKl(ll)/- 

6  DAK  1  (12)/  4.2B804170B26B9E-13/,  DAK  1(13)/ 

7  DAK  1  ( 1 4) /-l . 39140135641  18E- 14/ 

DATA  D  A  K  2  (  1 )  /  2 . 9333234388323E-0 1  /  , 

1  D A K 2  <  2 )  / -8 . 06 1 96784 7431 1 E -03/ ,  D A K 2 ( 3 )  / 

2  DAK.2  (  4  i  / -6 . 3229754885024E-04/ ,  DAK2<5)  / 

3  D A K 2 ( 6  i  / - 4 . 97457447684 06E-05/,  DAK2(7>  / 

4  DAK 2 ( 8 )  / -3. 495282 4044494E-06/,  DAK2(9)  / 

5  D A K 2  (  1 0 ) / - 2 . 4 4 7 3 2 6 7  1 5 2 1 3 7 E - 0 7 / ,  DAK2  ill)/ 

6  DAK2( 12) /-1 . 7271762150 154E-08/ ,  DAK2  (13)/ 

7  DAK  2 ( 14) / -  1 . 2324905529  1 55E -09/  ,  DAK2 (15)/ 

8  DAK2( 16) '-8. 392520 997724 0E-11/,  DAK2 (171/ 

9  DAK2( 18) /-6. 48013921 15345E-12/ ,  DAK2 ( 1 9> / 

1  DAK  2 ( 20 ) / -4 . 76  30382983364E - 1 3/ ,  D A K  2  <  2 1 ) / 


b (9 )  /-: 

b  ( 1 1 ; 


9  5  4  2  4  7  0  2  5 1  3  0  8  E  -  0  3  /  , 
431376  3  783919E-04/  , 
1631658269521E-06/  , 
7024377835588E-08/, 
0471633058877E-09/  , 
25358789454 13E- 12  /  , 
966296141  1602E-14/  , 
35798252031 10E-1 6/ 

6582062370270E-01 /  , 
634762881 0664E-02/ , 
4474440958980 E- 03/  , 
1066207768030E-05/ , 

1  5  2 1  8  8  0  7  2  7  7  1  6  E  -  0  7  /  , 
44284 02756595E-09/  , 
2080 2286 50 5 28E-  1 1  /  , 
6602263254  704E-13  /  . 
3616131657  0  0  2  E -  1  5 / 
576 47 27794 620E -03/  , 
3  5  9  1 6  0  B  0  2  6  8  8  2  E  -  0  7  /  , 
8682476755998E-08/  , 
2866061 8434  33E- 10/  , 
3235765530068E- 1 2/ , 

2  0  3  4  0  8  2  8  0  4  9  7  2  E  -  1  3  / , 
55656454580 15E- 15/ , 

529156918B258E-03/  , 
7131784226356E-06  /  , 
00100301 1 8434E-08/ , 
55098136 10 122E-1 1/  , 
83007 18795420E  - 12  /  , 
1052917947674E— 13/ , 
882051 1 1 00994E- 1 5/ , 


8 . 4  9  8  4  5  3  0  0  9  B  9  2  9  E  -  0  3  /  , 
•2. 700164898294 3E- 04  /  , 
3 . 0  2  3  9  7  7 1 2  4  0  9  5 1  E  -  0  6  /  , 
5 . 3  6 1 9  4  2  8  9  3  3  2  8  3  E  - 1 0  /  , 
7. 430 23834629 07E-1  1/  , 
4. 2381092575454E-13/  , 


2.42540 1 723331 4E- 
1.85786427751 1SE- 
1.3209068123950E- 
9. 24362451 07834E- 
6.4930783764891E- 
4.6072576360466E- 
3. 3062040943810E- 
2. 397733 1937830E- 
1 . 7551 01 3202373E- 
1 . 29498241  10081E- 
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A22 


2  D  A  K  2 (22) / - 3 . 

3  D  A  K  2 (24) /-2. 
DATA  D  A  K  3  i  1 ) 

1  D A K  3  ( 2 >  /  2. 

2  DAO  (4)  /  1. 

3  D  A  K  3  ( 6  )  / 

4  D  A  K  3  (  8  )  /  5. 
3  DAX'3  (10)/  6. 

6  DAO  ( 12)  /  1 . 

7  DAO  (  14)  /  2. 
DATA  DAJP(l) 

1  DAJP (2)  /- 1. 

2  DAJP  ( 4 )  /  1. 

3  DAJP  ( 6 )  / -8 . 

4  DAJP  (  8 )  /  2. 

5  DAJP  C 10) /  —  3. 

6  DAJP ( 12) /  4. 

7  DAJP ( 14) /— 3. 

8  DAJP  (16) 7  2. 

9  DAJP  < 18 ) /- 1 . 
DATA  DAJN  ( 1 ) 

1  DAJN  ( 2 )  /  8. 

2  DAJN  <  4 )  / -8 . 

3  DA J  N  ( 6  )  /  9. 

4  DA JN  ( 8 )  / - 4 , 

5  DAJN (10)/  9. 

6  DA JN  ( 1 2 ) / -  1 . 

7  PAJN< 14) /  1 . 

8  DAJN(16)/-1. 

9  D  A  J  N (  18) /  7. 

DATA  DA(l)  / 

1  DA (3)  / 

2  DA  (5)  /- 

3  DA (7)  /- 

4  DA (9)  / 

5  D  A  ( 1  1 )  /  - 

6  DAT  13) / - 

7  DA (15)/ 

DATA  DB  ( 1  )  /- 

1  DB  (  3  )  /- 

2  DB  (  5  >  / 

3  DB ( 7  >  /- 

4  D  B  (  9  )  /  - 

5  DB ( 1 1 ) / 

6  DB  < 1 3) /- 

7  DB  ( 1 5 ! / - 


5267962221 043E-14/ 

627869 1 4  34229E  -  1  5/ 

/  2. 8467582881  1 3 5 E - 0  1  /  , 
530730726 1 908E-03/  ,  DAK3 ( 3 i  /■ 
8490728394634E-06/,  DAK3(5)  /- 
05925634457 15E-09/,  DAK3(7)  /■ 
563576888 3  134E-1 1  /  ,  DAO(9>  /- 
8857 4 353 78 4 44E- 13/,  DAO  (11)/- 
173747626 172  IE- 14/,  DAK3 (13)/- 
5537477309706E-16/ 

/  6. 5321913131  1 4 6 E - 0 2 /  , 


DAK2 (23) /  9. 62005 1 5 1 58592E -  1 5 /  , 


4.83481  13033798E-05/  . 
1 . 0  1  4  1  8  4  9  1  1  7  8  5  8  E  -  0  7  /  , 

5 . 8  5  3  2  5  2  9 1 4  0  0  3  8  E  - 1 0  / , 

5 . 9  0  8  8  9  0  9  4  7  7  9  5  0  E  -  1  2  /  , 
8. 6858825645219E— 1  4/ , 
1. 6852314651 092E-15/, 


2026293368B82E-01/,  D  A  J  P ( 3  >  / 
6794842923050E-02/ .  DAJP (5)  / 
45560 29509887E-04/  ,  DAJP (7)  / 
2582786 094548E-05/,  DAJP ( 9 )  / 
76343991 13692E-07/ ,  DAJP (11)/ 
29611 3320030 1 E -09 / ,  DAJP (13) / 
5724588 1 36 190E-1 1/,  DAJP (15)/ 

2  6  1  2  0  6  5  3  0  9  5  7  7  E -  1 3 / ,  DAJP ( 17) / 
1260437448512E— 15  /  ,  DAJP (19)/ 

/  1 . 0859453963297E-02/ , 
53313194857 09E- 02/,  DAJN<3)  / 
7842072529426E-02/,  DAJN ( 5 )  / 
41674 06050 3 24E-03/,  DAJ N ( 7 )  / 

1  1  1573431 5683E  —  04/ ,  DAJN ( 9 )  / 
876336822084 0E-06/,  DAJN (11)/ 
507985 0013147E- 07/,  DAJN (13)/ 
5959991 741922E— 09/ ,  DAJN (15)/- 
2 4 1 0 3 3 5 0 5 0 0 3 0 E  —  1 1 / ,  DAJN ( 17) / 
39400971  15574E— 14/ ,  DAJN  (19)/- 
4. 9162732110 46 0E- 01/,  DA ( 2 >  / 
8. 2314 0 762854 08E-05/,  DA ( 4 )  /■ 
6.  1315888 053463E- 08/,  DA ( 6 )  / 
1. 819597 1537212E -09/,  DA ( 8 )  /• 
4. 5372404 34 204 2E- 1 1 / ,  DA  (10)/- 
3. 240891  1 933032E- 13/  .  DA(  12)  / 
2.  40 765247974 06E-1 4/  ,  DA  ( 1  4  )  / 
8. 179 00 7864774 OE-16/ 

2.  77571 35694423E-01 /  ,  DB  ( 2  >  / 
8. 4232852219009E-05/,  DB ( 4 )  /■ 
3. 4238972 021762E-07/',  DB  ( 6 )  /- 
2. 3637783684458E-09/,  DB(8)  / 
40995691 65819E-1 2/ ,  DB ( 1 0  > /- 
64874015137 02E-13/,  DB  ( 1 2  > / - 
772534306781  1  E -  1  4 /  ,  DB  ( 1  4 )  / 

1 179333757953E-16/ 


9. 7B01023626382E-03  /  , 
-1 . 97 1 46 1401 82 13E— 03  /  , 
9. 4288962 070198E-05  /  , 
-2. 2906787 091599E- 06/, 
3.4566393355956E-08/, 
-3. 5867369121499E  —  10  /  , 
2. 72696091 0663 4E-12/, 
-1.5876320 52383 OE-14/, 
7. 31 32752951537E  —  1 7 / 


-1 


152770681  1306E-01  /  , 
5  3  2  5  1 9  0  6  9  7  6  0  5  E  -  0  2  /  , 
32187 02601900E- 03/ , 
01297326891 35E-04, , 
8731 29698 1 239E— 06/ , 
3268766952539E-08/  , 
0766592266838E-10/ , 
3963176533104E-12/, 
32B8747562750E-15/ 

1  1  16493042749E— 03  /  , 
6176977617214E  —  06/ , 
8729580465652E-08/, 
4475282664204E-10/ , 
9965506584722E- 1 2/ , 
6209895256874E- 1 3/ , 
6938481 128449E-16/ , 


4. 4421283341992E-03/, 
2. 580403 184 1 87 1 E  —  06/ , 
6. 242868947 0973E- 09/, 

3.  1 699  1  04265667E- 1 0/ , 
5. 1 867422  109358E-1 2/  , 

4.  -yo  1 905766087 IE-  14/ , 
5 . 5  5  9  5  0  6  1  0  4  4  2  6  6  E  -  1  5  /  , 


I  F ( X . L T . 0 . i  GO  TD  300 
IF  (C. GT. 5.  )  GO  TO  200 
I F  (  X . GT . 1.2)  GO  TO  150 
T=  (X  +  X-l . 2) »C0N4 
T  T  =  T  -t-  T 
J  =N  1 

F 1 =AK  1  (J) 

F  2 = 0 . 
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TEMP  1  =  F  1 

Ft  =  TT*Fl-F2+AI03(J) 

F2  =  TEMP  1 
205  CONTINUE 

RTRX=SQRT ( R  X ) 

EC  =  EXP ( -C  i 

AI=EC* (T*F 1-F2+AK3 ( 1 ) )/RTRX 

J  =N  1 D 

F1=DAK3(J) 

F2  =  0. 

DO  206  1=1, MID 
J  =  J  - 1 
TEMP  1 =F  1 

F 1=TT*F 1 -F2+DAK3 ( J) 

F2  =  TEMP  1 
06  CONTINUE 

DAI=-RTRX*EC*(T*Fl-F2+DAK3(l)i 

RETURN 

300  CONTINUE 

IF (C. GT. 5. )  GO  TO  350 
T  = . 4  *C— 1 . 

TT=T+T 
J  =N3 

F  1  =  AJP  (  J  ) 

E 1 =A JN ( J ) 

F2  =  0. 

E2  =  0 , 

DO  305  1=1, M3 
J=J-1 
TEMP  1  =F  1 
TEMP2=E 1 

F 1  =  TT*F 1 -F2  +  AJP  ( J  ) 

E1=TT*E1-E2+AJN(J) 

F2=TEMP1 
E2=TEMP2 
305  CONTINUE 

AI  =  (T»E1-E2  +  AJN ( 1 ) ) - X #  ( T*F  1 -F2+-AJP  < 1 )  ) 

J  =  N  3  D 

F 1  =  DAJP  ( J ) 

E 1  =  D  A  J  N  ( J ) 

F  2  =  0 . 

E  2  =  0  . 

DO  306  1=1, M3D 
J  =  J-i 
TEMP  1  =F  i 
TEMP2  =  E  1 

FI  =  TT*F 1 -F2+DAJP ( J ) 

E 1  =  TT*E  1  -E2+-DAJN  ( J  ) 

F2=TEMP 1 
E2=TEMP2 
306  CONTINUE 

DAI  =  X*X*(T*F1-F2  +  DAJP(1) >+(T*El-E2+DAJN(l)  ) 
RETURN 
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350  CONTINUE  AIRY274 

T  = 1 0 . / C -  1 .  AIRY275 

TT=T+T  AIRY276 

J  =  N  4  AIRY277 

F 1 =A ( J )  AIRY278 

E 1 =B ( J )  AIRY279 

F2=0.  AIRY280 

E2=0 .  AIRY281 

DO  310  1=1, M4  AIRY282 

J  = J -  1  AIRY283 

TEMP  1 =F  1  AIRY284 

TEMP2  =  E  1  AIRY285 

F 1 =TT*F 1 -F2+A ( J )  AIRY286 

E 1 =  T T * E 1 -E2  +  B  ( J )  A1RY287 

F2=TEMP 1  AIRY288 

E2=TEMP2  AIRY289 

310  CONTINUE  AIRY290 

TEMP1=T*F1-F2+A(1)  AIRY291 

TEMP2  =  T*E1-E2  +  B(  1 )  AIRY292 

RTRX=SQRT (RX )  AIRY293 

CV=C-FPI  12  AIRY294 

CCV=COS(CV)  AIRY295 

SCV=S I N ( CV )  AIRY296 

AI= < TEMPI *CCV-TEMP2*SCV) /RTRX  AIRY297 

J  =N4D  AIRY298 

F 1 =DA ( J )  AIRY299 

El  =  DB  <  J )  AIRY300 

F2=0.  AIRY301 

E2=0 .  AIRY302 

DO  311  1=1, M4D  AIRY303 

J=J-1  AIRY304 

TEMPI =F 1  AIRY305 

TEMP2=E  1  AIRY306 

F 1 =TT*F 1 -F2+DA  ( J )  AIRY307 

E 1  =  TT*E 1 -E2+DB  < J )  AIRY308 

F2=TEMP1  AIRY309 

E2=TEMP2  AIRY310 

311  CONTINUE  AIRY311 

TEMP1=T*F1-F2+DA(1)  AIRY312 

TEMP2=T*El-E2+DB(l)  AIRY313 

E1=CCV»C0N5+.5*SCV  AIRY314 

E2=SCV*C0N5-.5*CCV  AIRY315 

DAI=(TEMPI*E1-TEMP2*E2)*RTRX  AIRY316 

RETURN  AIRY317 

END  AIRY318 

FUNCTION  6AMLN ( X )  GLN1 

c - -  -gi_N2 

C  WRITTEN  BY  D.  E.  AMOS,  SEPTEMBER,  1977.  SLN3 

C  REFERENCES  SLN4 

C  *  SAND-77-1518  GLN5 

C  *  COMPUTER  APPROXIMATIONS  BY  J.F.HART,  ET.Al.,  SIAM  SERIES  IN  6LN6 

C  APPLIED  MATHEMATICS,  WILEY,  196B,  P.135-136.  GLN7 

C  *  NBS  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS,  AMS  55,  BY  GLN8 

C  M.  ABRAMOWITZ  AND  I. A.  STEGUN,  DECEMBER.  1955,  P.257.  GLN9 

C  ABSTRACT  GLN10 
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GAMLN  COMPUTES  THE  NATURAL  LOG  OF  THE  GAMMA  FUNCTION  FOR  GLN11 

X.GT.O.  A  RATIONAL  CHEBYSHEV  APPROXIMATION  IS  USED  ON  GLN12 

8.LT.X.LT. 1000.  ,  THE  ASYMPTOTIC  EXPANSION  FOR  X.GE.1000.  AND  GLN13 
A  RATIONAL  CHEBYSHEV  APPROXIMATION  ON  2.LT.X.LT.3.  FOR  GLN14 

0.LT.X.LT.8.  AND  X  NON- 1 NTEGRAL ,  FORWARD  OR  BACKWARD  GLN15 

RECURSION  FILLS  IN  THE  INTERVALS  0.LT.X.LT.2  AND  GLN i 6 

3.LT. X.LT.8.  FOR  X  =  l. ,2 . 100. ,  GAMLN  IS  SET  TO  GLN  1 7 

NATURAL  LOGS  OF  FACTORIALS.  GLN 1 8 

GLN  1 9 

DESCRIPTION  OF  ARGUMENTS  GLN20 

INPUT  GLN21 

X  -  X.GT.O  GLM22 

OUTPUT  GLN23 

GAMLN  -  NATURAL  LOG  OF  THE  GAMMA  FUNCTION  AT  X  GLN24 

-  GLN25 


DIMENSION  GLN(IOO)  ,P(5) ,Q(2) , PCOE  ( 9 ) , QCOE ( 4 )  GLN26 

DATA  XLIM1 ,XLIM2,RTWPIL/  8.  ,  1000.  ,  9. 189385332047E-01 /  GLN27 

DATA  P ( l ) / 7 . 663451880000E-04/ ,  P  (2) /-5 . 940956 1 05200E-04/ ,  GLN28 

1  P (3) /7. 936431  104845E-04/  ,  P (4) /-2. 777777756577E-03/ ,  GLN29 

2  P(5)/8.333333333332E-02/  GLN30 

DATA  Q(l) /-2.777777777778E-03/,  Q  (2) /8. 333333333333E-02/  GLN31 

DATA  PCOE(l) /2.973786644B10E-03/,PC0E(2> /9.238194559028E-03/,  GLN32 

1  PCOE ( 3) / 1. 093 1 1595671 0E-01/,PC0E  (4) /3. 98067 1 3 1 0204E-0 1 / ,  GLN33 

2  PCOE (5 ) / 2 . 1 59943128461 E+00/,PCQE<6)/6. 3380679993B7E+00/ ,  GLN34 

3  PCOE (7) /2. 078247253179E  +  0 1 /  ,PCOE  (8) /3. 603677253002E+01 / ,  6LN35 

4  PCOE  <9 )/ 6. 2003838007 1 3E+01 /  GLN36 

GLN37 

DATA  QCOE ( l ) / 1 . OOOOOOOOOOOOE+OO/ , QCOE (2I/-8.906016659498E+00/,  GLN38 

1  QCOE (3 ) /9. 8225211047 14E+00/, QCOE (4) /6. 2003838007 13E+01 /  GLN39 

GLN40 

DATA  GLN ( 1 )  /O. 0/ ,  GLN (2)  / 0 . 0 / ,  GLN(3)  / 6 . 93 1 4 7 1 80559 9 E -0 1 / ,  GLN41 

1  GLN ( 4 )  /1.791759469228E+00/,  GLN (5)  /3. 1 78053330348E+00/ ,  GLN42 

2  GLN ( 6 )  /4. 787491 742782E  +  00/,  GLN  (7 )  /6. 5 7 9 25 1 2 1 2 0 1 0 E  +  0 0 / ,  GLN43 

3  GLN (8)  /8 . 525 161 36 1 065E+00  /  ,  GLN ( 9 )  /  1 . 060460290274E+0 1 / ,  GLN44 

4  GLN ( 1 0) / 1 . 2801 82748008E+0 1  /  ,  GLN ( 1  1  ) / 1 . 51 044 1 257308E+01 / ,  6LN45 

5  GLN ( 12) / 1 . 7502307845S7E+01 / ,  GLN ( 1 3 ) /  1 . 99872 1 449566E  +  0 1 / ,  GLN46 

6  GLN ( 14) /2. 255216385312E+01  / ,  GLN ( 1 5) /2. 51 9 1 221  1 8274E  +  0 1 / ,  GLN47 

7  GLN(lo)/2.789927138384E+01/,  GLN ( 1 7 ) /3 . 0671860 10608E+0 1 / ,  GLN48 

8  GLN ( 1 8 ) /3. 3505073450 1 4E  +  0 1  /  ,  GLN  ( 1 9 ) /3 . 639544520803E  +  0 1 / ,  GLN49 

9  GLN(20)/3. 93398841872 OE +01/,  GLN ( 2 1 ) /4 . 23356 1646075E+01 / ,  GLN50 

1  GLN ( 22 ) /4.538013889848E+01/,  GLN (23) /4 . 847 1 181 35 1 84E+0 l / ,  GLN51 

2  GLN (24) /5. 160667556776E+01/,  GLN ( 25 ) /5 . 4784729398 l 1 E+0 1 / ,  GLN52 

3  GLN  ( 26 )  /'5. 800360522299E  +  0 1  /  ,  GLN  ( 27 ) /6.  1 261701 76100E+0 1  /  ,  GLN53 

4  GLN (28) /6. 45575386270 1E+01/,  GLN (29) /6. 78B97431 3718E+01 / ,  GLN54 

5  G L N ( 3 0 ) / 7 . 1 2  5  7  0  3  8  9  6  7 1 7  E  +  0 1  /  ,  GLN (31 ) /7. 465823634883E+01 / ,  GLN55 

6  GLN  <  32) /7. 3 09222355 332E  + 01/,  GLN ( 33 ) /8 . 1 557959456 1 2E+0 1 / ,  GLN56 

7  GLN(34)/8. 50 544670 1758E+01/,  GLN(35) / 8. 85808 275422 0E+01/,  GLN57 

8  GLN(36)/9. 2136175 60 369E+01/,  GLN (37 ) /9 . 57 196945421 4E  +  0 1 / ,  GLN58 

9  GLN ( 38 ) /9. 933061 245479E  +  01 /  ,  GLN  ( 39 )/ 1 . 02968 1 996 1 45E  +  02/  GLN59 

DATA  GLN(40) / 1.0 663176026 06E+02/,  GLN ( 4 1 )/ 1 . t 03206397 1 48E+02/ ,  GLN60 

1  G L N ( 4 2 ) / 1  .  1403421  17815E  +  02  /  ,  GLN ( 43) / 1 . 1777 1 8B 13997E  +  02/ ,  GLN61 

2  GLN ( 44) /1.215330815154E+02/,  GLN ( 45 ) / 1 . 2531 727 1 1494E+02/ ,  GLN62 

3  GLN(46) /  1.291239 3 36391E  + 02/,  GLN  (47) /  1 . 329525750356E  +  02/ ,  GLN63 

4  GLN (48) / 1 . 368027226373E+02/  ,  GLN (491/1. 4 0673923648 2E +02/,  GLN64 

5  GLN (50) /l. 445657439463 E+02/,  GLN (5 1 )/ 1 . 4847776695 18E+02/ ,  GLN65 
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6  GLN  <  52  )  /  1 . 524095925845E  +  02/’ , 

7  G  L  N  ( 5  4  >  /  1 . 6  0  3  3 1  1  2  8  2 1 6  6  E  +  O  2  /  , 

8  GLN (56) /  1 . 683274454484E+02/  , 

9  GLN ( 58 ) /l. 763958484 070E+02/, 

1  GLN (60) / 1 . 8453382886 14E+0 2/ 

DATA  GLN (61 ) / 1 . 88628 1734237E  +  02/  , 

1  GLN (63) / 1 . 96866 1 81 6729E+02/ . 

2  GLN (65) /2.05168199  4826E  +  02/  , 

3  GLN ( 67 )  17. 1 353224 1 4946E+02/ , 

4  GLN (69) /2. 219564418191E+02  /  , 

5  GLN ( 7 i ) /2.304390435658E+02/, 

6  GLN (73)/ 2. 38978389561 8E+02/, 

7  GLN175) / 2. 47572914 0962E+ 02/, 

8  GLN (77) /2. 56221 1355500E+02/, 

9  GLN (79) /2.649216497986E+02/, 

1  GLN (8 1 ) /2. 736731 242857E+02/, 

2  GLN ( 83 ) /2.324742926876E+02/, 

3  GLN (85) /2.913239500943E+02/, 

4  GLN (87) /3.002209486470E+02/, 

5  GLN (89) /3. 09164 1935802E  +  02/ , 

6  GLN ( 9 1 ) /3. 181526396202E+02/, 

7  GLN (93) / 3. 27 1 852877038E  +  02/ , 

8  GLN(95)/3.362611819792E+02/, 

9  GLN(97)/ 3, 453794 070623E+02/, 
DATA  GLN (99) /3. 545390855 194E+02/, 

C  - - - 

5  NDX  =  X 

T=X-FL0AT (NDX ) 

IF (T.EQ.0.0)  GO  TO  51 
DX=XLIM1-X 


GLN ( 5  3 i /l. 5636083630 31E  + 02/ ,  GLN66 

GLN  (55) / 1.64320 1  1 22632E  +  02/,  GLN67 

GLN  (57) /I . 72352797 1392E+02/ ,  GLN68 

GLN  ( 59) / 1.804562914  175E  +  02/  ,  GLN69 

GLN  70 

G L N  ( 6 2 ) /  1 . 9 2 7 3 9 0 4 7 2 8 7 8 E  +  0 2 /  ,  GLN71 

GLN(64) /2. 010093163993E+02/ ,  GLN72 

GLN (66) /2.093425867525E+02/ ,  GLN73 

GLN (68) 12.  177369341  140E+02/ ,  GLN74 

GLN  (70)  /2. 2 6  1  9 0 5 4 8 3 2 3 7 E 0 2 /  ,  GLN75 
GLN (72) /2.347017234428E  +  02/  ,  GLN76 

G  L  N ( 7  4 ) /  2 . 4  3  2  6  8  8  4  9  0  0  3  0  E  +  0  2 / ,  GLN77 

GLN (76) /2.51B904022097E+02/ ,  GLN78 

GLN (78) /2. 6056494097 19E +02/,  GLN79 

GLN (80) / 2. 69291 09765 10E+02/  ,  GLN80 

GLN ( 82 )/ 2. 780675734404 E+02/  ,  GLN81 

SLN(B4> /2.368931332954E+02/  ,  GLN82 

GLN  ( 86 )/ 2. 9576660 13508E  +  02/ ,  GLN83 

GLN (88) /3.046868567657E+02/,  GLN84 

GLN  (90) / 3.  1365282994 99E  +  02/  ,  GLN85 

GLN (92) /3. 22663499 1267E+02/,  GLN86 

GLN (94) /3.317178871969E+02/,  GLN87 

GLN (96) /3.408150588708E+02/ ,  GLN88 

GLN (98) /3. 499541  1 8  0  4  0  8  E  +  0  2 /  GLNB9 

GLN  (1 00/3.59134  205  3696E  +  02/  GLN90 

GLN91 

GLN92 

GLN93 

GLN94 

GLN95 


IF1DX.LT. 0.0)  GO  TO  40  GLN96 

C  GLN97 

C  RATIONAL  CHE6YSHEV  APPROXIMATION  ON  2.LT.X.LT.3  FOR  GAMMA  ( X )  GLN98 

C  GLN99 

NXM=NDX-2  GLN 1 00 

PX=PCOE ( 1 )  GLN101 

DO  10  1=2,9  GLN 1 02 

10  P  X  =  T*P  X  +PCQE  ( I )  GLN 1 03 

QX  =  QCQE ( 1 )  GLN 1 04 

DO  15  1=2,4  GLN105 

15  QX  =  T*QX+QCOE  ( I )  GLN106 

DGAM=PX/QX  GLN107 

IF(NXM.GT.O)  GO  TO  22  GLN108 

IF (NXM.EQ.O)  GO  TO  25  GLN109 

C  GLN110 

C  BACKWARD  RECURSION  FOR  0.LT.X.LT.2  GLN 1 1  1 

C  GLN 1 1 2 

DGAM  =  DGAM/ ( 1 . +T )  GLN 1 1 3 

IF(NXM.EQ.-l)  GO  TO  25  GLN 1  14 

DGAM  =  DGAM/T  GLN  1 1  5 

GAMLN  =  ALOG ( DGAM )  GLN  1  16 

RETURN  GLN  1 1  7 

C  GLN  1  18 

C  FORWARD  RECURSION  FOR  3.LT.X.LT.8  GLN  1 19 

C  GLN 120 
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22  X X  =  2 . +  T 

DO  24  1  =  1 . N X M 
QGAM=DGAM*XX 

24  XX=XX+1. 

25  GAMLN=ALOG (DGAM) 

RETURN 

X.GT. XLIM1 

40  R  X  - 1 .  /  X 
RXX=RX*RX 

I F  < (X-XLIM2) . LT. 0. )  GO  TO  41 
P  X  =  Q ( 1 ) *RX  X  +  Q  <  2 ) 

GAMLN=F'X*RX+(X-.5)*AL0G(X)-X+RTWPIL 

RETURN 

X.LT. XLIM2 

41  PX=P ( 1 ) 

SUM=U-.5)*AL0G(X)-X 
DO  42  1=2,5 
PX  =  PX*RXX+P  ( I ) 

42  CONTINUE 

GAMLN=PX*RX+SUM+RTWPIL 

RETURN 

TABLE  LOOK  UP  FOR  INTEGER  ARGUMENTS  LESS  THAN  OR  EQUAL  100. 

51  IF  (NDX. GT. 100)  GO  TO  40 
GAMLN=GLN (NDX ) 

RETURN 

END 
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APPENDIX  B:  NOTATION 


Incident  wave  amplitude 

-l8Va 

Gravitational  acceleration 
Water  depth 

vcr 

Bessel  function  of  the  first  kind 

Wave  number 

Non-negative  integers 

Polar  coordinate 

Temporal  coordinate 

Horizontal  coordinate 

Horizontal  coordinate 

Bessel  function  of  the  second  kind 

Vertical  coordinate 

Incident  wave  angle 

Free  surface  displacement 

Polar  coordinate 

Angle  related  to  wedge  angle 

Value  related  to  wedge  angle 

Velocity  potential  function 

Horizontal  component  of  the  velocity  potential  function 
Defined  in  Equation  12 

Incident  wave  velocity  potential  function 
Reflected  wave  velocity  potential  function 
Scattered  wave  velocity  potential  function 

Finite  cosine  transform  of  <|> 

Wave  radian  frequency 


